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METHOD AND SYSTEM FOR OPTIMIZATION OF GENERAL SYMBOLICALLY 
EXPRESSED PROBLEMS, FOR CONTINUOUS REPAIR OF STATE FUNCTIONS, 
INCLUDING STATE FUNCTIONS DERIVED FROM SOLUTIONS TO 
5 COMPUTATIONAL OPTIMIZATION, FOR GENERALIZED CONTROL OF 

COMPUTATIONAL PROCESSES, AND FOR HIERARCHICAL META-CONTROL 
AND CONSTRUCTION OF COMPUTATIONAL PROCESSES 

CROSS-REFERENCE TO RELATED APPLICATION 
10 This application claims the benefit of provisional patent application nos. 

60/420,921 and 60/420,922, both filed October 23, 2002 

TECHNICAL FIELD 

The present invention is related to general optimization methods employed 

15 to determine relatively optimal solutions to various types of problems and to control of 
computational processes and construction of computational approaches to solving 
complex problems and, in particular, a general computational method and system for 
minimizing a general class of objective functions, for adjusting solutions to optimization 
over time, and for constructing and achieving relatively optimal computational control 

20 over a wide range of computational processes. 

BACKGROUND OF THE INVENTION 

Significant research-and-development effort is currently being applied to 
the mathematical and computational field of optimization. Optimization techniques are 

25 applied to a wide variety of everyday and complex theoretical problems in order to find 
solutions to those problems, including cost-effective and time-efficient control strategies 
for complex systems and organizations. For example, optimization techniques may be 
applied to schedule labor and equipment resources within a large manufacturing 
organization in order to most cost-effectively produce a given amount of manufactured 

30 product to meet a projected demand for the product at various points in time. As another 



2 



example, optimization techniques may be applied to a large control problem, such as the 
control of traffic lights within a city, in order to allow for as large a volume of traffic to 
move through the city as possible without generating unreasonable delays at various 
intersections and for those traveling on various secondary and tertiary routes within the 
5 city. Many natural systems can be viewed as seeking relatively optimal solutions to 
complexly defined problems governed by physical characteristics and principles and 
constrained by a large number of physical constraints. For example, complex chemical 
reactions can be considered to be governed by rather simple, thermodynamic principles 
by which individual components, such as molecules, seek minimum energy states and 

10 maximum entropy, constrained by various ways in which molecules can reconfigure 
themselves and exchange energy amongst themselves according to quantum mechanics. 
Even more complex biological systems can also be considered to be governed by 
chemical and thermodynamic principles as well as by analogous principles involving 
information content and exchange, compartmentalization and modular design, and other 

15 considerations. Thus, optimization problems may encompass an extremely wide variety 
of mathematically and computationally expressed models of natural phenomena, design, 
organization, and control of complex systems, and organization, transmission, and 
processing of information. 

Current approaches for finding near optimal and optimal solutions for 

20 mathematically modeled problems are limited. When the number of decision variables 
and constraints in such problems increases from the small number of variables and 
constraints normally employed in simple, textbook problems to the large number of 
variables and constraints normally encountered in real-world systems, the computational 
resources required for seeking near optimal and optimal solutions increases exponentially 

25 in most cases. Current techniques cannot be satisfactorily applied to any but the simplest 
types of problems. Many currently available techniques involve applying oppressive 
constraints to optimization models, such as requiring variables to be continuous and 
requiring the hyper-dimensional volume representing the set of possible solutions of 
optimizations problems, or problem domain, to be convex. 
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Researchers, developers, system modelers, and investigators of many 
different types of complex system behavior have therefore recognized the need for more 
generally applicable and efficient methods and systems for optimization. 

5 SUMMARY OF THE INVENTION 

Various embodiments of the present invention include methods and 
systems for generic optimization of problems by an approach to minimizing functions 
over high-dimensional domains that mathematically model problems for which near 
optimal or optimal solutions are sought. These embodiments receive a mathematical 

10 description of a system, in symbolic form, that includes decision variables of various 
types, including real-number-valued, integer-valued, and Boolean-valued decision 
variables. The objective function with high dimensional domain that is minimized in 
order to find near optimal and optimal solutions may be accompanied by a variety of 
constraints on the values of the decision variables, including inequality and equality 

15 constraints. Various embodiments of the present invention incorporate the objective 
function and constraints into a global objective function. The local and global minima of 
this function occur at critical points. One or more of these critical points is a near optimal 
or optimal solution for the modeled problem. Various embodiments of the present 
invention transform the global objective function and a procedure for finding critical 

20 points into a system of differential equations in terms of continuous variables and 
parameters, so that powerful, polynomial-time methods for solving differential equations 
can be applied for identifying critical points of the function. These embodiments 
undertake an interior-point-based method and employ a global gradient-descent field 
formulation in order to efficiently traverse the hyperdimensional domain of the global 

25 objective function to find local minima. Near a local critical point, where the global 
descent field is ineffective in providing better approximants of the variables, 
embodiments of the present invention provides a steering method that generates a local 
descent-field across the local critical point towards the solution. Once the vicinity of the 
local critical point is left behind, the global gradient descent method is resumed, and 

30 various embodiments of the present invention employ a method to establish a reasonable 
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direction for resuming a descent path through the hyper-dimensional domain toward local 
minima. One of these local minima represents near optimal solutions of the problem. 
Embodiments of the present invention also provides for distribution and decomposition 
of the global gradient descent-field and local gradient descent-field optimization methods 
5 using multiple threads and agents, respectively, in order to allow for parallel computation 
and increased time efficiency. Various embodiments of the present invention further 
include approaches for adjusting solutions to optimization problems relatively 
continuously in time, without needing to recompute the optimization solution de novo. 
While many embodiments of the present invention are specifically directed to various 
10 classes of optimization problems, other embodiments of the present invention provide a 
more general approach for constructing complex hierarchical computational processes 
and for optimally or near optimally controlling general computational processes. 

BRIEF DESCRIPTION OF THE DRAWINGS 
15 Figures 1A-B illustrate a very simple, two-dimensional optimization 

problem. 

Figure 2 is a control-flow diagram for a highest level view of an 
optimization method and system that represents one embodiment of the present invention. 

Figure 3 illustrates a hypothetical problem used to illustrate modeling and 

20 preprocessing. 

Figures 4-7 illustrate further aspects and characteristics of the hypothetical 
problem using the illustration conventions of Figure 3. 

Figure 5 illustrates costs associated with shipping goods from a 
manufacturing plant to warehouses and from warehouses to customers. 
25 Figure 6 illustrates additional characteristics of the hypothetical problem. 

Figure 7 shows a simple two-manufacturing-plant, three-potential- 
warehouse, five-customer hypothetical problem associated with numerical costs, 
capacities, and demands. 

Figure 8 is a more detailed, control-flow diagram for the preprocess step of 
30 the high-level control-flow diagram shown in Figure 2. 
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Figure 9 shows an initial symbolic, mathematically expressed model for 
the specific optimization problem illustrated in Figure 7. 

Figure 10 shows transformation of the hypothetical mathematical model 
into a model employing only float variables. 
5 Figure 11 shows the mathematical model for the hypothetical problem 

transformed to standard form. 

Figures 12A-D illustrate the effect of binding variables on transformed 
equality constraints. 

Figure 13 shows the symbolic, mathematically expressed model for the 
10 hypothetical problem following transformation of equality constraints to inequality 
constraints. 

Figure 14 graphically illustrates the effect of the binding variables. 

Figure 15 shows the mathematical model for the hypothetical problem 
following introduction of the envelope variable z. 
1 5 Figures 1 6A-D illustrate the concept of barrier functions. 

Figure 17 shows the symbolic, mathematical model for the hypothetical 
problem following addition of barrier constraints. 

Figure 18 illustrates the use of a gradient field within a problem domain in 
order to provide a direction for an optimization trajectory leading from an initial point to 
20 a local minimum. 

Figure 19 illustrates an optimization trajectory. 

Figure 20 graphically illustrates an iterative optimization method based on 
a near-optimal state-vector trajectory y and near-optimal control policies w*and r . 

Figure 21 provides additional, graphical illustration of the overall iterative 
25 optimization method. 

Figure 22 is a control-flow diagram that describes the iterative 
optimization step (203 in Figure 2). 

Figure 23 shows a surface boundary for a bowl-shaped problem domain 
with a number of prominent central protrusions. 
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Figure 24 shows the same problem domain as shown in Figure 23 with 
critical points explicitly indicated. 

Figure 25 shows a magnified view of a critical point. 

Figure 26 illustrates a first-search technique. 
5 Figure 27 illustrates a second strategy for searching a problem domain. 

Figure 28 illustrates the modules and modular organization of an 
optimization thread. 

Figure 29 illustrates, in part, the contents of the database maintained by the 

super module. 

10 Figures 30-34 illustrate, using the illustration conventions of Figure 28, the 

flow of information among the hierarchically organized modules of an optimization 
thread during one outer-loop iteration of the optimization method. 

Figure 32 illustrates data flow during a next step of the inner loop (2205 in 

Figure 22). 

15 Figure 33 indicates that the date-exchange steps illustrated in Figures 31 

and 32 are successively repeated in a number of inner-loop iterations. 

Figure 34 illustrates the data exchange involved in costate vector 
calculations undertaken by the p module 2810 in step 2209 of Figure 22. 

Figure 35 graphically illustrates computation of the h> vector. 
20 Figure 36 illustrates multiplication of the Hessian matrix by the w vector 

to produce the G vector in a distributed fashion. 

Figure 37 shows a decomposed multiplication of the Hessian matrix by the 

w vector. 

Figure 38 shows a mathematical manipulation of the computation of a G- 
25 vector partition that facilitates decomposition of the Hessian-matrix inversion by y- 
module agents. 

Figure 39 illustrates computation of a H>-vector partition w 1 by the second 
>>-module agent. 

Figure 40 is a control-flow diagram illustrating a highest conceptual level 
30 of the super module. 
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Figure 41 is a control-flow diagram for the routine "initial point" called in 
step 4004 of Figure 40. 

Figure 42 is a control-flow diagram for the routine "choose_next_T." 

Figure 43 is a control-flow diagram of the routine "evaluate," called in 
5 step 4012 of Figure 40. 

Figure 44 is a control-flow diagram for the routine "report." 

Figure 45 is a control-flow diagram for the routine 
"integrate_w_and_y_to_T 5 " called in step 4010 of Figure 40. 

Figure 46 is a control-flow diagram for the routine "integrate y." 
1 0 Figure 47 is a control-flow diagram for the routine "compute Ax," called in 

step 4625 of Figure 46. 

Figure 48 shows that the Hessian matrix can be decomposed by spectral 
decomposition. 

Figure 49 shows successive state-function calculations by an optimization 
1 5 method for a control problem. 

Figure 50 illustrates a shortened interval between calculations of the state 
functions illustrated in Figure 49. 

Figure 51 illustrates an analogy for consequence arising from repeated, 
current-state-based optimizations. 
20 Figure 52 illustrates the concept of a sliding window. 

Figure 53 illustrates a state function calculated in a prior window extends 
from time Ti-AT to time TV AT. 

Figure 54 illustrates extension of the previously calculated state function 
to extend over the time interval of the current window. 
25 Figure 55 illustrates consideration, at a current point in time, of past events 

and environmental changes that may result in a previously calculated state function being 
out-of-date, and not optimal, at the current point in time, as well as an approach to 
incorporating information from the previous state function in finding a new, optimal state 
function for the current frame of reference represented by the current window. 
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Figure 56 shows repair of the previously calculated state function in order 
to provide a near optimal state function for a current time window. 

Figure 57 illustrates a single level of meta-control of a computational 

process. 

5 Figure 58 illustrates multiple levels of meta-control and/or hierarchical 

construction of a complex computational process. 

DETAILED DESCRIPTION OF THE INVENTION 

The present invention relates to the generic optimization of 

10 mathematically-formulated problems, relatively continuous adjustment of solutions to 
optimization control problems, and to a more abstract, higher-level hierarchical 
construction of, an optimal or near-optimal control of computational processes. Various 
embodiments of the present invention include complex software programs running on a 
single-processor computer system, or running, in parallel, on multi-processor computer 

15 systems, or on a large number of distributed, interconnected single-and/or-multiple 
processor computer systems. These programs are quite mathematical in nature, and 
involve elaborate and difficult-to-describe mathematical techniques. For this reason, the 
present invention is described in part, below, with reference to a concrete problem, with 
reference to many hundreds of equations, and with reference to many graphical 

20 illustrations and control-flow diagrams. Although mathematical expressions, alone, may 
be sufficient to fully describe and characterize the programs that implement processes that 
represent various embodiments of the present invention to those skilled in the art of 
software development, the more graphical, concrete-problem-oriented, and control-flow- 
diagram approaches included in the following discussion are intended to illustrate various 

25 embodiments of the present invention in a variety of different ways that may be 
particularly accessible to readers with particular backgrounds. The present invention is 
described below in the following subsections: (1) Optimization; (2) Repair; and (3) 
Meta-Control and Hierarchical Process Construction. 
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Optimization 

Figures 1A-B illustrate a very simple, two-dimensional optimization 
problem. Figures 1 A-B show two different views of a two-dimensional surface, with the 
surface mathematically expressed as follows: 
5 z = f (xl,x2) 

where the domain of the problem is the two-dimensional surface and the goal of the 
problem is to identify values of xl and x2 for which the dependent variable z is minimal. 
In other words, the graph of the surface is a collection of points {xl, x2, z) characterized 
by a function of two variables, xl and x2, plotted with respect to the x/-axis 102, the x2- 

10 axis 104, and the z-axis 106. As a concrete example, z might represent the cost associated 
with manufacturing widgets of length xl and width x2. An optimization problem with 
regard to this model may be to find the optimal widget length and widget width in order 
to reduce the cost of manufacturing a widget to the lowest possible value. In most cases, 
a closed-form expression for the solution, or even a mathematical model, for an 

15 optimization problem is not available. Instead, one may begin with some initial, feasible 
values for the decision variables, in the present case xl and x2 values corresponding to a 
point on the two-dimensional surface 110, and then attempt to find a path, denoted in 
Figure 1A by the dark curve 112, from the initial, feasible position 110, to a local or 
global minimum value of the function z min that occurs at a point 1 14 on the surface with 

20 coordinates (xl min , x2 min , and z min ). 

Note that one possible approach for finding a near-optimal or optimal 
solution to the problem would be to compute z for every possible pair (xl, x2). However, 
in general, even for extremely simple problems, such an approach would be prohibitively 
computationally. In many cases, the domains for decision variables, such as decision 

25 variables xl and x2, are subsets of the real numbers, or subsets of other natural domains, 
which may greatly decrease the potential number of decision-variable-values for which 
the function would need to be evaluated in order to conduct a complete, combinatory 
search. However, even in these cases, the combinatory search method is generally not 
feasible. 
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It should be noted that, although many two-dimensional-surface 
illustrations, such as Figures 1 A-B are used in the following discussion, they are provided 
only for illustration purposes. In real-world problems, there may be thousands, hundreds 
of thousands, or millions of variables and constraints, leading to extremely high- 
5 dimensional problem domains. In general, these problem domains may be considered to 
be hyper-dimensional volumes, or manifolds, bounded by hyper-dimensional surfaces 
specified by the constraints. It is not possible to illustrate such problem domains, but 
techniques used to address hyper-dimensional problem domains may analogized to 3- 
dimensional illustrations. For example, a feasible point in a two-dimensional problem 

10 domain is a point of a two-dimensional surface, while a feasible point in a hyper- 
dimensional problem domain is a point within the hyper-dimensional volume, or 
manifold, enclosed by the hyper-dimensional boundary surface. 

Many other approaches are possible for traversing the domain of a 
mathematical model from an initial, feasible point to a local or global minimum point. 

1 5 One can, for example, compute the gradient for the function at each point, and proceed 
along the surface in the direction of the gradient in order to most steeply descend along 
the surface to a local or global minimum. This approach is best suited for convex 
problems, in which a single global minimum exists. The hyper-dimensional volume 
domain for a real-word optimization problem is quite often non-convex, and includes 

20 many critical points, including saddle points, local maxima, and local minima in addition 
to a global minimum. Therefore, a steepest-descent approach using gradients may be 
used to traverse the hyper-dimensional problem domain to a local minimum, but a 
searching method must nonetheless be employed in order to visit and evaluate multiple 
local minima. For example, in Figure 1A, had the trajectory dictated by the steepest- 

25 descent method ended at the local minimum represented by point 116, then some other 
approach would be needed in order to detect and evaluate the global minimum at point 
112. Figure IB shows the problem-domain surface and path to the global minimum 
shown in Figure 1 A rotated by approximately forty-five degrees about the z-axis. 

Figure 2 is a control-flow diagram for a highest level view of an 

30 optimization method and system that represents one embodiment of the present invention. 
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In a first step, 201, a particular system or problem to be optimized is mathematically 
modeled, and the mathematical model serves as input to subsequent optimization steps. 
Many of the embodiments of the present invention include graphical user interfaces, text- 
based interpreters, and other such interfaces in order to facilitate mathematical modeling 
5 of problems. Below, a particular mathematical model for a hypothetical problem is 
provided as an example of a mathematical model for an optimization problem. In 
general, the mathematical model is provided in symbolic form, comprising a set of 
mathematical expressions including symbolic terms and operators. Next, in step 202, the 
method and system undertake a preprocessing step in order to transform the symbolically 

10 expressed mathematical model to a standard form needed by methods that represent 
embodiments of the present invention. Preprocessing steps are discussed below with 
respect to the hypothetical mathematical model, and the preprocessing steps are formally 
described. Next, in step 203, the symbolic standard form of the mathematical model is 
used to set up an iteration and to iteratively optimize the function obtained using an 

15 instance of the interior-point method and a gradient-steered descent-field trajectory in 
order to identify one or more local minima within a hyper-dimensional problem domain. 
Various strategies of descent field mechanisms are employed to find and evaluate a 
possibly large number of local minima in order to identify a near-optimal solution to the 
mathematically modeled problem. Various embodiments of the present invention employ 

20 different search strategies. It is important to note that neither the methods that represent 
the embodiments of the present invention, nor any other optimization methods currently 
employed, are guaranteed to always find a global minimum within a given tolerance or a 
reasonable amount of time using a reasonable amount of computer resources. However, 
most practical optimization problems are bounded, in that the number of local minima is 

25 finite. Finally, once a best solution is found within some specified tolerance or the 
maximum amount of time using a specified amount of computing resources, this solution 
is transformed, in a post-processing step 204, in order to provide one or more near- 
optimal solutions in the context of the mathematical model for the problem originally 
provided in step 201. 
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The mathematical modeling step and preprocess step, 201 and 202 in 
Figure 2 5 respectively, are perhaps best appreciated using a concrete, hypothetical 
problem. Figure 3 shows a hypothetical problem used to illustrate the mathematical 
modeling and preprocessing steps. The hypothetical problem involves p manufacturing 
5 plants, 302-304, w warehouses, 306-311, and c customers, 314-318. A manufacturing 
plant, such as manufacturing plant 302, can ship goods to any of the w warehouses, 
indicated in Figure 3 by various types of lines, such as dashed line 320 connecting 
manufacturing plant 302 and warehouse 306. Each warehouse, in turn, can ship goods to 
any of the c customers, with paths between warehouses and customers indicated by 

10 various types of lines, such as the dash line 322 interconnecting warehouse 306 and 
customer 314. In Figure 3, a different style of line is used to show the path from a 
particular manufacturing plant to the various warehouses, and from a particular 
warehouse to the various customers. 

Figures 4-7 illustrate further aspects and characteristics of the hypothetical 

15 problem using the illustration conventions of Figure 3. For example, the 
notation x£ w indicates the amount of goods shipped from manufacturing plant i to 

warehouse j and terms of the form xj£ indicate the amount of goods shipped from 

warehouse j to customer k. Although there are w potential warehouses, not all w 
warehouses need to be used. A Boolean variable b j indicates whether or not warehouse j 

20 is employed, with the Boolean value "1" or "TRUE" indicating that the warehouse is 
used, and the Boolean value "0" of "FALSE" indicating that the warehouse is not used. 
In Figure 4, for example, warehouses 1 (306), 3 (308), and 5 (310) are used, as indicated 
by the solid borders of the squares representing used warehouses and the Boolean values 
"1" associated with the used warehouses, and warehouses 2 (307), 4 (309), and all 

25 warehouses following warehouse 5, including warehouse w 311, are not used. The used 
warehouses are alternatively numbered, in Figure 4, with sequential numbers in 
parentheses so that the three used warehouses are sequentially numbered 1, 2, and n, 
where n is the total number of used warehouses and equals 3 in the example shown in 
Figure 4. 



13 



Figure 5 illustrates costs associated with shipping goods from 
manufacturing plants to warehouses and from warehouses to customers. A term of the 
form indicates the per item cost of shipping goods from a manufacturing plant i to 

warehouse j\ and a term of the form cj£ indicates the per item cost of shipping goods 

from warehouse j to customer k. In Figure 5, each of the shipping routes is labeled with 
the corresponding per-item cost. 

Figure 6 illustrates additional characteristics of the hypothetical problem. 
With each manufacturing plant, there is an associated manufacturing capacity, with a term 
having the form j. indicating the shipping capacity of manufacturing plant /. Each 
warehouse has a receiving capacity, expressed by a term v y indicating the capacity of 

warehouse j to receive goods. Also, there is a cost of operation for each warehouse, with 
a term in the form r y indicating the cost of operating warehouse j. Finally, there is a 

demand of goods associated with each customer. A term of the form ^indicates the 
demand for goods by customer k. In Figure 6, the supply capacities, receiving capacities, 
cost of operations, and demands associated with manufacturing plants, warehouses, and 
customers are all explicitly detailed. Finally, Figure 7 shows a simple two- 
manufacturing-plant, three-potential-warehouse, five-customer hypothetical problem 
associated with numerical costs, capacities, and demands. 

The optimization problem illustrated in Figures 3-7 is to figure out the 
lowest cost configuration of operating warehouses, manufacturing-plant-to-warehouse 
shipping volumes, and warehouse-to-customer shipping volumes such that all the 
customer demands are met and no capacity is exceeded. A more formal, symbolic 
mathematical expression of the optimization problem follows. 

First, the total cost associated with any particular system configuration is 

expressed as: 

total a*- zl>r*f + ZI>; c *; c+ 2>a 

,=1 y=! >=1 *=1 7=1 
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In other words, the total cost is the sum of the costs associated with all shipments from 
manufacturing plants to warehouses, from warehouses to customers, and the costs of 
operating warehouses. A minimum is sought in the hypothetical optimization problem 
for the total cost. In addition, the model includes various constraints. First, a 
5 manufacturing plant cannot ship a greater amount of goods than the capacity of the 
manufacturing plant, expressed as follows: 

7=1 

Customers all receive at least an amount of goods equal to their demands, expressed as 
follows: 

7=1 

Furthermore, the amount of goods sent by a manufacturing plant to each warehouse is 
less than or equal to the warehouse's capacity for receiving goods, expressed as follows: 

Similarly, a warehouse ships an amount of goods less than or equal to the warehouse's 
15 capacity to receive goods, expressed as follows: 

k=l 

Warehouses do not store goods. Each warehouse ships out everything that it receives, 
expressed by the conservation constraint that follows: 

20 The shipment pathways illustrated in Figures 3-7 are unidirectional or, in other words, 
goods are not returned, expressed as follows: 
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xf >0 and xj k c >0 

For the hypothetical problem it is assumed that there is sufficient manufacturing capacity 
and warehouse-receiving capacity to satisfy the customer demands, expressed as follows: 

p w c 

i=\ 7=1 *=1 

Goods may be shipped to and from only operating warehouses, expressed by the 
following two constraints: 

*f(l-&,) = 0 

The number of operating warehouses is less than or equal to the total potential number of 
warehouses w: 

Finally, the terms yj are Boolean values, expressed as follows: 

be {0,1} 

The above mathematical model may be slightly transformed to more 
concisely express the cost function and associated constraints. First, the amount-of- 
goods-shipped parameters of the types x^andxjmay be replaced by scaled shipment 
terms as follows: 

m P» = _v_ 9 o</wf <1 



X 

<;=-^-,o<<<i 
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The warehouse-capacity-based constraints, provided above, turn out to be redundant, and 
can therefore be removed from the set of constraints without affecting the problem model 
or solutions obtained for the model. The conservation constraint indicating the 
warehouses do not store goods is modified to use scaled shipment terms, as follows: 

±smf >£d k m] k < 

1=1 *=1 

Non-operable warehouses do not ship goods, as expressed by the following constraint: 

mj<bj 

The total cost may be expressed, using the scaled shipment of goods terms, as follows: 

/=! y=l y=l *=1 j=\ 

The capacity constraints are expressed, in terms of the scaled shipment of goods terms, as 
follows: 

The above-listed capacity constraints can be recast as an inequality with respect to 0, as 
follows: 

± Si rf -±d k m% >0 

/=! k=\ 

The fact that the number of operating warehouses must be less than or equal to the total 
potential number of warehouses, that goods cannot be returned, and that the bj terms are 
Booleans, are now represented as: 
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,/»;;> o 

6,e{0,l} 

Figure 8 is a more detailed control-flow diagram for the preprocess routine 
invoked in step 202 of the high-level control-flow diagram shown in Figure 2. Each step 
5 of the preprocess routine illustrated in Figure 8 is further discussed below with reference 
to corresponding transformations of the hypothetical mathematical model described 
above. First, in step 801, the preprocess routine receives the initial, symbolic, 
mathematically expressed model for the optimization problem, including constraints. In 
step 803, the preprocessing routine attempts to remove as many redundant constraints as 

10 possible. An example of the removal of redundant constraints occurs in the above 
transformed model, with the warehouse-capacity constraints in the original form of the 
model recognized as redundant and removed from the transformed model. Next, in step 
805, all Boolean and integer-value decision variables, such as the Boolean variables bfs, 
in the hypothetical example, are transformed into floating-point variables. In step 807, 

15 the resulting symbolic, mathematically expressed model is placed in a standard form. In 
step 809, equality constraints are replaced by relaxed inequality constraints. Each equality 
constraint is relaxed with three inequality constraints using an additional positively 
valued variable. In step 81 1, an additional variable z is added, along with an additional 
constraint to convert the criterion into an additional constraint. Finally, in step 813, the 

20 transformed optimization problem is converted into an Interior Point Method barrier 
function by adding together terms. Each term is an evaluation of a barrier function in 
terms of an expression representing a constraint multiplied by a control variable. This 
Interior Point Method Barrier function contains two sets of control variables u and r. As 
a result of the preprocessing routine, the mathematical model is symbolically expressed in 

25 a form that can be optimized by the iterative optimization technique shown as step 203 in 
Figure 2. 

Figure 9 shows the initial symbolic, mathematically expressed model for 
the specific optimization problem illustrated in Figure 7. The cost function, 902, is to be 
optimized with respect to vectors of integer variables x and Boolean variables b. The 
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constraints include manufacturing-plant-capacity constraints 904, customer-demand 
constraints 906, conservation constraints 908, constraints 910 that model the fact that no 
flow can occur to or from a warehouse that is not selected by the optimization, that is the 
corresponding Boolean variable is zero, a constraint that only two of the three warehouses 
may be operational 912, a constraint specifying that there are no return shipments of 
goods from customers to warehouses or from warehouses to manufacturing plants 914, 
and the constraint that the b/s are Boolean-valued variables 916. 

The third step of the preprocessing routine, step 805 in Figure 8, involves 
transforming Boolean, integer, and other non-floating-point variables to floating-point 
variables. Each Boolean variable 6, is replaced with a floating-point value y t in the range 
[0,1] and with an additional constraint: 

y t <y-y i ) = o 

The additional constraint implies that the floating-point variable >>, can only have one of 
the two values 0 or 1. Each integer variable is replaced by an equivalent base-2 
polynomial with k terms, where k is approximately equal to the base-2 logarithm of the 
maximum value that the integer variable can have. Each term of the base-2 polynomial 
comprises a power of 2 multiplied by a Boolean coefficient, and thus the base-2 
polynomial is, in turn, replaced with k floating-point variables y 0 , yu... yk-h all within the 
range [0,1] along with k constraints in the form: 

The conversion of decision variables of various types to floating-point numbers is 
summarized formally as follows: 



19 



x i g R=> x, is float 

x, ; <£ R — > x ; is transformed to floats: 

Vx, e Bool, x, => y t e [0, 1], y, (1 - y t ) = 0 

Vx # e Z + , x,. => J] 6,. 2' = b 0 2° + 6, 2 1 + b 2 2 2 + . . . + b k _ x 2 h ~ l 

where maxvalue(x l ) = 2* - 1 
log 2 (maxvalue(x,)) = A: 

^0^1^2-^-1 G [0,1] 



Figure 10 shows transformation of the hypothetical mathematical model into a model 
employing only float variables. 

Next, the mathematical model is placed in standard form. The formal 
definition for a standard-form symbolic problem representation is provided as follows: 

minimize C(x) 
subject to : 

g / (x)>0/or/ = 0roM 1 -l 
h. (x) = 0 for i = M x toM 2 -\ 



where x = 



and x 0 e ffi 



In other words, a standard form model includes a cost function C(jc), a number Mi of 
inequality constraints in the form 

&00*o 
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and a number of equality constraints of the form: 

jj(x)=0 

Figure 1 1 shows the mathematical model for the hypothetical problem transformed to 
standard form. Note that transformation of the original equality and inequality constraints 
to constraints with 0 on the right-hand side, as needed for the standard form, is a simple 
matter of adding terms to both sides of the constraint equations that are not already in 
standard form. 

The next step in the preprocessing routine (809 in Figure 8) is to change 
all equality constraints to inequality constraints. In order to do so, each equality 
constraint is transformed to two inequality constraints, each inequality constraint 
involving the left-hand expression of the original equality constraint and a positive 
binding variable. Transformation of an equality constraint h/x) to two inequality 
constraints is formally expressed as follows: 



*.(x) 



/*, (x) + s y >0 
-h J (x) + s J >0 



15 In addition, the constraint 

Sj >0 

is added to ensure that the constraints expressing h/x) are feasible. 

Thus, following transformation of equality constraints to inequality 
constraints, the total number of constraints increases by a number 2k equal to the original 

20 number of equality constraints in the standard-form model of the problem. Figures 12A- 
D illustrate the effect of binding variables on transformed equality constraints. Figure 
12A shows the original inequality constraint h/x) with the value 0 (1202). When the 
equality constraint h/x) is transformed to three inequality constraints, as discussed above, 
then the range of values for h/x) expands about 0 to include all values between -sj 1204 

25 and Sj 1206. Thus, transforming an equality into three inequality constraints essentially 
relaxes the equality constraint from a point to an interval. Relaxing equality constraints 
facilitates the iterative optimization methods, to be described below, by preventing 
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equality constraints from overconstraining, or dominating, the initial descent towards one 

or more local minima. Relaxing equality constraints also facilitates an interior-point 

descent, since equality constraints force optimization trajectories to problem-domain 

boundaries. However, as the optimization begins to converge, the relaxation of equality 

5 constraints needs to be removed. This is accomplished, as shown in Figures 12C and 

1 2D, by decreasing the absolute value of the barrier variable Sj in order to more and more 

tightly constrain the value of h/x) to a small region symmetrically disposed about 0. The 

degree of relaxation of an equality constraint can be expressed by a coefficient r, for each 

binding variable sj added to the criterion. Transformation of equality constraints to 

10 inequality constraints can be formally expressed as: 

minimize C(x) 
constrained by : 

g,(x)>0 / = 0,l,...,m-l 

hj(x) = 0 j = w,l,...,/w + k-l 

u 

2k-\ 

minimize C(x) = C(x) + ]>] r i s i 

/=o 

constrained by : 

g,(x)>0 i' = 0,l,...,m + 2*-l 

Figure 13 shows the symbolic, mathematically expressed model for the hypothetical 
problem following transformation of equality constraints to inequality constraints. Figure 
14 graphically illustrates the effect of the binding variables r h In Figure 14, the 

15 optimization method begins at initial point 1402 and descends through a problem-domain 
1404 bounded by the cup-shaped surface 1406 towards a minimum 1408. Initially, the 
coefficients r, for the binding variables s t are small, providing a large total range of 
values, shown in Figure 14 by the disk-shaped region 1410, about the initial point 1402. 
As the descent through the problem domain approaches the minimum 1408, the 

20 coefficients r, of the binding variables s f are increased, more and more tightly 
constraining the solution sets until, near the optimum 1408, the transformed equality 
constraints again become extremely tightly constrained 1412 because as the optimization 
method proceeds the s t becomes smaller and smaller. 
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The next step of the preprocessing routine involves adding an envelope 
variable z along with an additional constraint to the mathematical model of the 
problem (811 in Figure 8). The envelope variable z is defined as follows: 

z>C(x) 

5 This leads to the following inequality constraint: 

z-C(x)>0 

Following introduction of the envelope variable z, the formal statement of the 
problem becomes: 

minimize z 

with respect to x 0 , ., x n _ l9 s 0 , s, ,...,s k _ ]9 z 
constrained by : 

g y ( x ) ' = Oto/w-l 

gj(x,Sj) y =0to2£-l 

z-c(x)>0 

10 Figure 15 shows the mathematical model for the hypothetical problem following 
introduction of the envelope variable z. 

Figures 16A-D illustrate the concept of barrier functions. Assume, for 
example, that the problem domain for an optimization problem is the wedge-shaped 
volume with a parabolic base 1602 shown on Figure 16A. An initial point for an 

15 optimization 1604 has a relatively large z- value 1601, which decreases to 0 at the 
minimum 1608. In Figure 16B shows an x-y plane projection 1610 of the problem 
domain 1602. The optimization trajectory from the initial point 1604 to the minimal 
point 1608 is shown as three line segments 1612-1614. In an interior point method, 
when traversing the problem domain, the optimization trajectory must remain within 

20 the problem domain. In other words, the current point in a descent trajectory through 
the problem domain to a minimum must always remain within the problem domain, 
interior to the boundary of the manifold in a hyper-dimensional problem domain. For 
example, as shown in Figure 16B, the initial trajectory 1612 is angled towards the 
edge of the volume corresponding to the problem domain, and would intersect the 

25 boundary at point 1616 and, if continued in the same direction, would result in 
passing through the boundary to points exterior to the problem domain along the 
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dashed line segment 1618. In order to keep the trajectory of the optimization within 
the problem domain, it is necessary to include barrier functions within a problem 
model that have little or no effect within the problem domain, but rise to large values 
as the optimization trajectory nears the boundary surface of the problem domain, in 
5 order to deflect the trajectory back from the boundary into the problem domain. For 
example, in Figure 16B, as the optimization trajectory 1612 approaches pointl616, 
the barrier functions added to the mathematical model begin to rapidly increase in 
order to drive the value of the envelope variable z to high values and therefore deflect 
the optimization trajectory away from the inner surface of the problem domain. 

10 Following deflection, a new optimization trajectory 1613 is established which leads, 
through continued optimization, to a new second approach to a boundary point 1620, 
just prior to which the barrier constraints again force the value of the objective 
function z to large values in order to again deflect the optimization trajectory back 
into the interior of the problem domain and to a new, final trajectory 1614. 

15 Figures 16C-D illustrate the operation of the barrier functions. 

Assuming a current optimization trajectory 1622, as shown in Figure 16C, consider 
the value of the objective function z with respect to the x-coordinate, shown in Figure 
16D, for a traversal of the problem of the domain along the current trajectory 1622 
shown in Figure 16C. As can be seen in Figure 16D, the value of z is relatively low 

20 for values of x within the interior of the problem domain along the current trajectory 
1622. However, as x decreases towards the left boundary point 1624 or increases 
towards the right boundary point 1626, the barrier functions have the effect of rapidly 
increasing the value of z 1628 and 1639, respectively. In other words, for each 
inequality in the model a term involving a barrier function is added to the interior 

25 point objective function in order to constrain trajectories to the interior of the problem 
domain. If a trajectory leaves the problem domain, or reaches the exterior of the 
enclosing surface of the problem domain, the optimization methods that represent 
some of the embodiments of the present invention could fail. 

Adding terms involving barrier functions to the model is a simple way 

30 to ensure that the optimization-trajectory stays within the problem domain. A barrier 
function is added for each constraint, since it is the constraints that specify the surface 
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enclosing the problem domain. A suitable barrier function term for the optimization 
techniques employed in various embodiments of the present invention is the negative 
of a positive variable coefficient w, the control variable associated with the term times 
the logarithm of the constraint expression g„ or, in other words: 
5 -w,log(g,) 

As the value of the constraint expression g, decreases towards 0, the value of the 
logarithm of the constraint log(g,) decreases rapidly with a corresponding rapid 
increase in absolute value. By taking the negative of the log (g,), a rapidly increasing 
value is obtained as the constraint expression g, closely approaches 0, approaching 

10 infinity in the limit, thus preventing the constraint from ever decreasing to 0. The 
variable coefficients w, can be considered to be control variables, much like the 
coefficients r, of the above-described binding variables. The w, and r, coefficients are 
modified, during an optimization, in order to steer the optimization trajectory in a 
most efficient manner towards local minima. However, unlike the n coefficients, the 

15 Ui coefficients are initially large, and are decreased by one of the embodiments of this 
invention, as the optimization trajectory approaches a local minimum, since local 
minima often fall close to, or on the surface of, the problem domain. A formal 
expression the addition of barrier functions to the objective function is provided as 
follows: 

i=m+2k 

20 minimizeF(x,s,z,u,r) = z- ^ M /l°g(&) 

Figure 17 shows the symbolic, mathematical model for the hypothetical problem 
following addition of barrier constraints. The addition of the barrier functions 
completes the preprocessing steps embodied in the preprocessing routine in a first, 
described embodiment. 

25 Next, an overview of the core optimization method that represents one 

embodiment of the present invention is provided. Figure 18 illustrates the use of a 
global gradient descent field within a problem domain in order to provide a direction 
for an optimization trajectory leading from an initial point to a local minimum. At 
each point within the problem domain, such as point 1802, the partial derivative 

30 vector of the objective function F(x,s,z,u,r) may be calculated as a vector of first 
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partial derivatives of the objective function with respect to the decision variables jc„ 
the binding variables s i9 and the envelope variable z. Similarly the Hessian of the 
function F(x,s,z,u,r) given as the matrix of second partial derivatives of 
F(x,s,z,u,r) with respect to the decision variables jc„ the binding variables s„ and 
5 the envelope variable z can be computed at every point of the domain The gradient of 
the objective function is defined as the vector at each point in the domain computed 
as the product of the inverse of the Hessian and the partial derivative of.F(x,s,z,u,r) . 
It provides a direction 1804 from the point 1802 representing a descent direction in 
the value of the objective function with respect to the decision variables and binding 

10 variables. Since a gradient can be calculated for each point within the problem 
domain, the gradient of the objective function specifies a vector field, pointing in a 
direction of descent and, specifying a, path to a minimum within the problem domain. 
In general, gradients for points, such as point 1802 in Figure 18, far from local 
minima have relatively large absolute values, while gradients for points close to a 

15 local minima, such as point 1806, have gradients with relatively small absolute 
values. 

The mathematical description of the optimization methods is 
facilitated by a number of notational conventions. First, the objective function 

F(x,s,z,u,r) 

20 can be expressed as: 

F(y,u,r) 

with 

x 



= y 
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U m+2k J 



= U 




Thus, the objective function can be thought of as a function of a vector y of decision 
5 variables, binding variables, and the envelope variable, the vector u 9 a vector of 
barrier-function coefficients, and the vector r, a vector of binding variable 
coefficients. The vector y is referred to as the "state vector," and the vectors u and r 
are referred to as "control vectors." The vector y has N elements, where N = n 
(decision variables) + k (binding variables) + 1 (envelope variable). 

10 Figure 19 illustrates an optimization trajectory. The trajectory 1902 is 

plotted with respect to the lengths of vectors r and u as horizontal axes, and the value 
of the objective function F as a vertical axis. The optimization trajectory begins at an 
initial point 1904 and descends in the value of the objective function F towards a 
local minimum. During the descent, the value of the barrier functions expressed by 

15 the length of vector u, decreases towards 0, while the value of the coefficients of the 
binding variables, as expressed as the length of vector r, increases. The trajectory is 
divided into segments of length Ax, such as a first segment 1906. During 
optimization, the vectors y and u are continuously recalculated to drive the value of 
the objective function lower. In Figure 19, the value of the binding variable 

20 coefficients, represented as vector r, is also shown to be recalculated during each Ax 
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interval, although in the described embodiment of the present invention, the binding 
variable coefficients are adjusted only towards the end of an optimization trajectory. 

The discrete optimization technique may be expressed as follows: 

y* + i=y*+^(y*>r„ii,) 

5 where k = 0, 1, kconvergence- In other words, a new state-vector value yk+i is 
computed, in each iteration of a computational optimization method, from a current 
state-vector value and a discrete function <p k that depends on j;*, r ki and u k , with the 
iteration count k ranging from 0 to a value k convergence at which point the state vector 
for a local minimum of the objective function might be determined. In a first step of 

10 optimization, the state and control vectors, y, «, and r may be parameterized by a 
continuous iteration parameter x, where x ranges from 0 to an iteration horizon T. By 
parameterizing vectors y, w, and r by the continuous iteration variable x, the iteration 
above is transformed to a first order differential equation in which the y variable is 
the dependent variable and «, and r are control variables. Well-known differential- 

15 equation-based mathematical technique can be used for computing efficient 
optimization trajectories as solutions of this differential equation. This is discussed 
below. The objective function is now expressed formally as: 

Hy(*)'"{*)>r(*)) 

where 0 < r < T . For example if iteration y M = y k +<f> k (y k ,r k9 u k ) is the Newton 
20 iteration, the corresponding differential equation has the following form 



dr 



It should be noted that Newton's method is but one approach to 

construct the differential equation providing an expression for — . It was chosen for 

dr 
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simplicity and for ease of computation. However, any of many other approaches may 

ds 

be used in embodiments of the present invention for approximating — . 

dx 

In the following discussion, the partial derivative vector G and Hessian 
matrix H for the objective function F are frequently used. The partial derivative 
5 vector G is defined as follows: 



GMr).«(r),r(r)) = 



dy(r) 



dF(y(r),u(T),r(r)) 
dF(y(T),u(T),r(r)) 

dF{y{r),u{r),r{r)) 

Qy N -i (^) 



10 The Hessian matrix H is a matrix of size AfrAf, defined as follows: 
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ac 0 


5G, 


dG 2 


3G„_, 












5G, 


5G 2 


dG^., 




3* 


3* 





SG 0 9G, 3G 2 dCJ w _, 
d>v, 5y w _, 



dy 0 
dG, 

9G-, 



9G, 



gG 0 

aG, 

5G, 



3G, 



5y 0 ^1 



ac 0 

5G, 
5G 2 



^N-l 
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d 2 F(y(T),u(T),r(T)) 



d 2 F(y(T),u(T),r(r)) 



d 2 F(y(T),u{T),r(r)) 



dy N _ x dy N _, 



The normalized rate of change of the state vector y with respect to the continuous 
iteration variable x, 



when expressed in terms of the Hessian matrix H and vector G, or as: 

dy = w 
dr p 

where w = H" l G and p = G r H -1 G . 

The optimization problem can be recast into a minimum time problem 
which is a special case of a generic, optimal control problem. In other words, by 
continual ization of the objective function F (y, u, r) , continuous methods can be used 

to obtain a computational approach for computing an optimization trajectory y{r) 
consisting of state- variable changes that bring the value of the objective function 
F(y,u,r)from an initial point to a minimum. However, the problem domain is 

enormous, and hyper-dimensional, and so, without further computational constraints, 
or meta-level constraints, although a method for descending the gradient field within 
the problem domain has been obtained, the problem of computationally descending 
the gradient field is slow convergent. 



dy 



H G 



G r H'G 
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An optimal control problem is formulated. It seeks to speed-up the 
convergence by minimizing the convergence horizon and satisfying the requirements 
on the control variables. In order to control minimization of the objective function 
F(^(r),w(r),r(r)) the optimal control formulation is used to steer, using the 

control variables u(r),r(r), the optimization trajectory by selecting values for the 

control variables that minimize the convergence horizon. The optimal control 
problem for controlling the optimization trajectory is formally expressed as follows: 

minimize j{u{r) 9 r{z) 9 T) = pr + ^G r QG 

with respect to u(r),r(r) } re[0,r], where, Q is a symmetric positive definite 
matrix, subject to the constraints: 

(1) forr e [0,r] 

dr p 

(2) «,e(0, Mmax ) 

(3) «,.(r + Ar)<u,(r) 

(4) n e^.oo) 

(5) /;(r + AT)ar ( (r) 

Minimization of the objective functional for the minimum convergence horizon 
problem, J"(u(r),r(r),r), forces F(y,u,r) to converge to a minimum, where G = 
0, in the smallest computational horizon T, where the convergence interval is the 
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integral of the continuous iteration variable. The solution to the optimal control 
problem provides an optimal control policy w*(r) and re[0,r] that 

generates a state-vector-change optimization trajectory y (r) . The first term in the 

integral defining the objective function J r (u(r) 5 r(r),r) forces the optimization 

5 trajectory to an efficient as possible trajectory y* with respect to the continuous 
iteration variable r . The second term within the integral forces the vector G towards 
0, thereby forcing the objective function F to a minimum. The combination of these 
terms therefore produces an optimal control policy u(x) and r(x) to force the 
objective function F to a local minimum through an optimal state- vector-change 

10 trajectory y*(r). Now the constraints of the optimal control problem (l)-(5) are 
discussed. The first constraint involves the normalized steepest descent for the state- 
vector-change trajectory, and is the primary vehicle by which the original 
optimization problem is incorporated into the new minimum-time control problem. 
The second and third constraints ensure that the control vector values u remains non- 

1 5 negative and less than some maximal control value u max and that the w- vector values 
do not increase during the optimization process. The fourth and fifth constraints are 
analogous to the second interior constraints with respect to the binding variable 
coefficients r,-. The fourth and fifth constraints ensure that the binding variable 
coefficients r, do not fall below a minimum value r mim and that the binding variable 

20 coefficients n do not decrease during the optimization process. 

A primary result of optimal control theory, known as the Pontryagin 
"Minimum Principle," provides a framework for the minimum convergence horizon 
problem to produce an optimal control policy u(x) and r(x) r e [0,r] for finding an 

optimal state-vector-change optimization trajectory y(x). This approach involves the 
25 Hamiltonian function, defined as follows: 



The Hamiltonian has a minimum value as a function of u and r at u(x) and r*(x). In 
other words, the Hamiltonian has a minimum value for the optimal control policy. 
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When the control policy w*(r) is optimal, there exists an absolutely continuous 
function p(r) , called the costate, such that the derivatives of the state and costate 
functions y(x) and p(r) are as follows for all r <= [0,T] : 



dy* _ dH(y*(r) 9 u(T) 9 p(T)) 
dr dp 

dp = r a^(y'(T) t n'(r) > p(r)) " ir 
dr dy 



5 satisfying the boundary conditions: 



y(o) = y. 

p (y) _ dQT (y( T ), uCT), r(T))QG(y(T), u(T), r(T)) 



forallu(r),W(y*(r),u , (r),p(r))<H(y*(r),u(r) ) p(T)). 
The control problem solution is summarized as follows: 



dy' _ w 
dr p 

dp* _ d ( w 

y(o) = y 0 

p(T) = H(T)QG(T) 

r w(y'(r),u(T),r'(T)) ^ T w(y (r), M (r),r(r)) 
P p(y(r),u'(r)y(r))- P p(y {r),u{r),r{r)) 
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The Hamiltonian is analogous to the Hamiltonian of a mechanical 



system, which expresses the total energy of a mechanical system. The function y*(r) 
is analogous to specifying the lowest energy positions of the components of a 
mechanical system in a continuous time interval, and the costate function p*(x) is 
5 analogous to specifying the lowest energy momentum of components of a mechanical 
system in a continuous time interval. Using the above-derived results allows for 
computation, in a forward direction, along the optimization trajectory, of the partial 
derivative 



can be computed backward from x = Tfrom the forward-computed y(x). Again, with 
analogy to a classical mechanical problem, the momentum of the components of a 
mechanical system are essentially determined from the change in position of the 
15 components, over time. Thus, momentum, or costate, needs to be computed 
backwards from forward computed positions. 



and produces the optimal control policy when the Hamiltonian is minimized. This 
minimization problem can also be solved using Newton's approximation method, as 



dr 



1 0 The partial derivative 



dp 



Minimization of the Hamiltonian is itself an optimization problem, 



20 follows: 



= 0 and 



en 



= o 



du (r) 



dr(r) 



35 











~d 2 ft 


d 2 n~ 




~dH 










du 2 


dudr 




du 






'to. 




d 2 H 


d 2 ft 




dn 










_drdu 


dr 2 J 




. dr 



(y(r k ) 9 u(T k ) 9 r(r k ) 9 p(r k )) 



where 



*(y( T t)>uM,r(T k ) 9 p(T k )) = 



d 2 n d 2 n 



du 2 

d 2 n 



dudr 

d 2 n 



drdu dr 2 



dft 

du 
dft 

dr J 



{yM>uM>r(r k ),p(T k )) 



Thus, formulation of the minimum time problem and optimal control problem 
provides a method for computing state vector y and control vectors u and r in a 
forward direction, and computing the costate vector p in a backward direction, over a 
range of the continuous iteration variable x between 0 and a horizon T. This then 

10 leads to an iterative optimization method by which the objective function F is forced 
to a local minimum via an optimal state- vector-trajectory y*(x) under an optimal 
control policy u *(r) and r*(x). 

Figure 20 graphically illustrates the iterative optimization method 
based on the above-derived state-vector-change optimization trajectory y* and optimal 

15 control policies u* and r*. The iterative optimization method can be thought of as 
comprising two nested loops, an outer loop controlled by the iteration variable K and 
the inner loop controlled by the iteration variable k. The outer nested loop represents 
integration of the partial differentials of the state-vector-trajectory and control policy 
over a continuous iteration range of 0 to T in order to find a state vector y that 

20 minimizes the objective function F. Within each such outer-loop iteration, multiple 
forward computations of the control policy vectors u and r and the state vector^ over 
small intervals of the continuous iteration variable r are carried out in steps r k , where 
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k goes from 0 to a final k value, or, in other words, T 0 ,T l9 x 2 ,... 9 x ftnai , and then, at 

x final > T , the results of the computations ofy(xk) and u(x0 are used to back calculate 

the costate vectors p(tt) from t = T to x = 0. During each inner loop, the computation 
of the y(x0 values are based on prior computation of u(xk) and r(xk), and the 
5 computation of both the y(x0 and jufx/J are carried out using costate values p(x0 
computed at the end of a previous iteration of the outer loop. In the first iteration, the 
control costate values are approximated as a constant value. 

In Figure 20, the outer loop is represented by the outer circle 2002, and 
the inner-loop computations are represented by step-like radial and circumferential 

10 segments, such as step 2004. The iterative optimization process begins with initial 
values for the state vector, control vectors, and costate vector of y 0 ,u 0 ,r 0 , and p 0 , 
respectively 2006. The outer-loop iteration variable K is set to 0 (2008). Then, a 
number of step-like, inner-loop iterations are carried out, beginning with iteration 
k = 0, each successive iteration of the inner loop incrementing the inner-loop variable 

15 A: by one. In the first, inner-loop iteration, control vectors u(xk+i) and rfa+i) are 
computed 2004. Then, the value of the state vector at x = x k+I , y(x k +i), is computed by 
a number of successive computations 2010-2016. At the end of the first iteration of 
the inner loop, a determination is made as to whether a local minimum has been 
reached. If so, then the iterative optimization method terminates, indicated in Figure 

20 20 by the dashed line 2018. Otherwise, a next inner-loop iteration is carried out 
2020-2021. Inner-loop iterations continue until x k > T 2022, at which point the 
costate vector values p(x) are computed in a number of steps 2024-2030 from x = Tto 
x = 0, and the outer-loop iteration control variable K is incremented to conclude a full 
outer-loop iteration and return to the beginning point 2032 for a next outer-loop 

25 iteration. 

Considering Figures 19 and 20 together provides a graphical, 
conceptual view of the iterative optimization process. Each inner-loop iteration of 
the optimization process can be considered to correspond to a segment in the 
optimization trajectory shown in Figure 19, such as segment 1906. One full outer- 
30 loop iteration, from x = 0 to x = T can be considered to produce a large segment of the 
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optimization trajectory, such as the optimization trajectory segment beginning at 
initial point 1904 and extending down to point 1908. Multiple outer-loop iterations 
therefore compute a number of large optimization trajectory segments that eventually 
drive the state vector to a value at which the objective function F is at least a local 
5 minimum. Figure 21 provides additional, graphical illustration of the overall iterative 
optimization method. As shown in Figure 21, a first outer-loop iteration produces the 
optimization trajectory segment from an initial point 2102 to point 2104. A second 
outer-loop iteration computes the optimization trajectory segment 2106 from point 
2104 to point 2108. Third outer-loop iteration produces the optimization trajectory 

10 segment 2110 from point 2108 to local optimum 2112. Note that the value of the 
horizon, T, may change with each outer-loop iteration. Thus, in Figure 21, the 
endpoint of each optimization trajectory segment is marked T K , with K = 1,2, and 3, 
indicating that the endpoint of each outer-loop iteration serves as the initial point for 
each subsequent outer-loop iteration It should also be noted that the magnitude of Ax 

1 5 for each inner-loop step may also vary. 

Figure 22 is a control-flow diagram that describes, in general, the 
iterative optimization step (203 in Figure 2). First, initial values for y, u, r, and p are 
chosen, the outer-loop iteration variable K is initially set to 0, and an initial horizon T 0 
is computed in step 2202. Next, in an inner loop comprising steps 2204-2208, the 

20 control vectors u and r and state vector y are propagated forward from x = 0 to x = T K 
in x increments of r k . First, in step 2204, a next value for the control vectors u and r 
is computed. Then, in step 2205, the state vector^ is propagated from x = Xk to x = 
Xk+Axk in step 2205. In step 2206, a determination is made as to whether the 
optimization has converged. If so, then the iterative optimization has completed, and 

25 returns. Otherwise, a determination is made as to whether Xk+Ax remains less than 
Tk, in step 2207. If not, then x k +j is assigned to x k +Ax, and the inner-loop control 
variable k is incremented, in step 2208. However, if x k is greater than or equal to T K , 
then the costate values are computed backwards from T K to 0, in step 2209, and the 
iteration variables and a database containing the computed values of the state vector, 

30 control vectors are updated in step 2210. 
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In convex problems, in which there is a single critical point, the 
iterative optimization method described above with reference to Figures 19-21 
normally can converge in a single thread within a reasonable amount of time even for 
extremely large hyper dimensional problem domains with thousands of decision 
5 variables. However, optimization problems in applications are usually not convex, 
and the problem domains may feature many critical points. Figure 23 illustrates one 
such problem domain. In Figure 23, the surface boundary for the problem domain 
2303 is bowl-shaped with a number of prominent central protrusions, such as 
protrusion 2304. These protrusions may be thought of as a small mountain range 

10 with three peaks (2304, 2306 and 2308) and a deep local valley, or bowl 2310. The 
global optimization amount for the problem domain lie deeply below the mountain 
range in the depths of the outer bowl. 

Figure 24 shows the same problem domain as shown in Figure 23, or 
critical points explicitly indicated. The peaks of the mountain range 2304, 2306 and 

15 2308 represent local maxima within the problem domain, and the low point in the 
small valley 2310 represents a local minimum. Saddle points within the mountain 
range, such as saddle points 2312 and 2314, also represent critical points. Figure 25 
shows a magnified view of a critical point, such as critical point 2304 in Figure 24. 
As shown in Figure 25, the partial differentials calculated for a number of points 

20 2502-2504 near the critical point in the jc and y directions calculated for a number of 
points 2502-2504 near the crucial point 2506 are all small in absolute value. 
Therefore, the gradient calculated for the points 2202-2304 near the critical point 
2506 are also close to 0. To recapitulate, the derivative of the objective function with 

' dF(y(T)Mr),r(T)) ) T 

25 optimization trajectory near a critical point. Therefore, when the gradient is 0 or near 
0, as it is for points near a critical point, then the rate of the state-vector trajectory 
with respect to x is also 0. Moreover, G is sensitive near the critical point, since small 
local variations can cause relatively large magnitude changes in G and can even 
switch the sign of G. However, when G is small, then Newton's method, and other 

30 similar methods used to compute state-vector trajectories, are generally not useful for 



respect to x, 



= G, gives a means for determining the 
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finding an optimization trajectory. In other words, the gradient field generally 
provides the direction for near-optimal or optimal trajectories, and therefore, while 
the gradient field is 0, all possible trajectories appear to be equally good. It is a 
situation analogous to attempting to use a compass at the magnetic north or south 
5 poles. As will be discussed later, a different method for computing the optimization 
trajectory for the state vector needs to be employed near critical points. 

There is another consideration with regard to critical points. This 
consideration involves search-like strategy that needs to be employed near critical 
points by a global optimization method in order to attempt to identify and evaluate a 

10 number of potential minima within a problem domain. Figure 26 illustrates a first- 
search technique. In Figure 26, an iteration of the outer-loop of the optimization 
method to horizon T n has carried the state vector to point 2602. This point is near the 
critical point 2306, the peak of the mountain range. The fact that the gradient is 
nearly 0 at this point informs the optimization method that it is arrived at the vicinity 

15 of a critical point. However, its next step depends on the type of critical point at 
which it has arrived. In the case of a local maximum, such as point 2306, the 
optimization method needs to employ some method for perturbing the optimization 
trajectory away from the local maximum 2306 to continue a descent 2308 towards a 
local minimum. However, in the neighborhood of a local minimum, such as local 

20 2310, the optimization method needs to employ a technique to ascend out of the local 
minimum 2312 to another point 2314 from which the optimization method can 
continue to descend 2316 towards a local minimum. 

Figure 27 illustrates a second strategy for searching a problem domain. 
In Figure 27, the search method has again arrived at a point 2602 near a local 

25 maximum 2306. However, upon detecting the local maximum, the optimization 
method spawns a number of additional, relatively independent optimization threads 
which continue descending in different directions 2702-2705. One of the relatively 
independent threads 2704 reaches a point 2708 near a saddle point 2314. At this 
point, this thread, 2704 forks, or spawns, an additional thread, so that each of the 

30 resulting relatively independent threads descends from the saddle point in different 
directions 2710 and 2712. Another of the originally spawned, relatively independent 
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threads 2702 reaches a local minimum 2310. At this point, rather than spawning new 
optimization threads, the relatively independent thread 2702 terminates, reporting the 
local minimum back to its parent thread 2701. Parent threads can select the lowest 
local minimum identified by their children and return that local minimum to their 
5 parents, resulting in a recursive search of the problem domain to identify a near 
optimal or optimal solution to the optimization problem. In alternative embodiments, 
parent threads can use information about reported local minima to truncate the 
problem domain, in order to modify the problem domain based on the local minima, 
and communicate the truncation of the problem domain to all other relatively 

1 0 independent threads. 

By either the single-threaded technique shown in Figure 26, or the 
multi-threaded technique shown in Figure 27, the optimization method needs to have 
a means of determining the type of critical point that is arrived at when the gradient 
falls to a value near 0. The critical-point beacon for illuminating types of critical 

15 points is the Hessian matrix, described above. When the gradient is near 0, then 
when all elements of the Hessian matrix are negative, a local maximum has been 
reached. When the gradient is 0 and all elements of the Hessian matrix are positive, 
then a local minimum has been reached. Otherwise, when the gradient is 0, and the 
elements of the Hessian matrix have mixed signs, then an inflection point or saddle 

20 point has been reached. The critical-point beacon properties of the non-singular 
Hessian matrix are formerly described as follows: 



for a point y , if G | y « 0 

if all eigenvalues of H are negative, then y is near a local or global maximum 
if all eigenvalues of Hare positive, then y is near a local or global minimum 
otherwise y is near a saddle point 



25 A single optimization thread, as discussed above with reference to 

Figures 26 and 27, comprises a collection of hierarchically organized modules that 
may execute concurrently in separate computational environments in order to 
distribute the large computational tasks involved in propagating the state vector 



41 



according to the above-described optimization method. Figure 28 illustrates the 
modules and modular organization of an optimization thread. The highest-level 
module is referred to, in this document, as a "super" 2802. The super includes a 
relatively large database 2804 of intermediate results. The super is responsible for 
5 computing each successive horizon T K , for choosing an initial point, in the case of the 
root thread in a multi-threaded optimization method, as described above with 
reference to Figure 27, and for launching and coordinating the activities of 
subordinate modules. The first-level subordinate modules include a ujmd r module 
2806, a y module 2808, and a p module 2810. The w and r module 2806 is 

10 responsible for the forward computation of the control vectors u and r, as shown in 
step 2204 in Figure 22. 

The y module 2808 is responsible for propagation of the state vector 
through successive intervals Ax, shown as step 2205 in Figure 22. The p module 
2810 is responsible for computing costate vector values backwards from T k to 0 at the 

15 end of each outer-loop iteration of the iterative optimization method, shown as step 
2209 in Figure 22. 

The w_and_r module 2806 distributes the control vector forward 
calculation amongst m agents 2812-2815. The y module 2808 and the p module 2810 
each distributes its respective task amongst n agents 2816-2818 and 2820-2823, 

20 respectively. Distribution of the y propagation by y module 2808 over the w^-module 
agents 2816-2818 is an important approach to decomposing the overall optimization 
computations and providing opportunities for parallel computing. The state-vector 
propagation computations undertaken by they module 2808 and its respective agents 
2816-2818 are the central, and most computationally expensive, tasks involved in 

25 computing the state-vector-change optimization trajectory. In an alternative 
embodiment of the present invention, the p module uses an efficient, interpolation 
method for estimating costate vector values, and thus does not need to distribute the 
problem among n p-module agents. Thus, the /^-module computes costate without 
distribution. In other alternative embodiments, the /^-module may distribute costate 

30 interpolation over a small number of p-module agents. 
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The database 2804 maintained by the super 2802 is large, and may be 
contained not only in memory, but in mass storage devises directly or remotely 
coupled to one or more computer systems over which the modules of a thread are 
distributed. Figure 29 illustrates, in part, the contents of the database maintained by 
5 the super in various embodiments of the present invention. In general, the super must 
maintain the values of various quantities calculated by the modules at each 
subinterval x k of the horizon from 0 to T. Thus, in Figure 29, each sub-column of the 
first column of data structures, such as sub-column 2902 of the vector G data 
structure 2904, represents the value of the corresponding computed entity at a 

10 particular value of x. In Figure 29, sub-column 2902 represents the value of the 
vector G at x = 0. The super stores values of the state vector 2906, the control vectors 
u and r 2908-2909, the vector h> 2910, the scaler p 291 1, a Morse-flags vector 2912, 
the costate vector p 2914, the scaler values of the objective function F 2915, the 
vector G 2904, the value of the Hamiltonian 2916, the current Tand iteration control 

15 variables and other such local control variables 2918, and the inverses of the Hessian 
matrix 2920-2927. 

Figures 30-34 illustrate, using the illustration conventions of Figure 
28, the flow of information among the hierarchically organized modules of an 
optimization thread during one outer-loop iteration of the optimization method. In 

20 the case of a root thread, as shown in Figure 30, an initial point for the optimization 
trajectory is determined 3002. In addition, initial values for the state vector, costate 
vector, control vector, and control variables, among other data structures, are 
determined for initializing the various data structures and elements 3004. Non-root 
threads inherit a starting point and current values for the various computational 

25 entities, such as the state and control vectors, from parent threads. As shown in 
Figure 31, in the first step of the inner loop (2204 in Figure 22), the w_and_r module 
2806 receives the values of the state vector, costate vector, vector m>, and the inverse 
of the Hessian from the super 3102 and distributes these values and/or portions of 
these values to wandr-module agents 3104-3107, receives intermediate computed 

30 results back from the «_and_r-module agents 3108, 3109, 31 10, and 3111, computes 
a new control vector value, and returns the new control vector value 3112 back to the 
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super. Figure 32 illustrates data flow during the next step of the inner loop (2205 in 
Figure 22). In this step, the y module 2808 receives the newly computed control 
vectors u and r from the super 3202, distributes these values, along with additional 
values obtained from the super and/or computed by the y module, including the 
5 Hessian, vector other such values, to the y agents 3204-3206, receives intermediate 
computed values back from the y agents 3208 3209 3210, and computes and returns 
to the super computed values for the state vector, w vector, differential of the w vector 
with respect to the state vector, and the inverse of the Hessian, among others 3210. 
The y module carries out a number of individual propagation steps for each inner- 

10 loop iteration, as discussed above. Figure 33 indicates that the data-exchange steps 
illustrated in Figures 31 and 32 are successively repeated in a number of inner-loop 
iterations. Figure 34 illustrates the data exchange involved in costate vector 
calculations undertaken by the p module 2810 in step 2209 of Figure 22. The p 
module receives computed values of the state vector, control vectors, and other such 

15 values from the super 3402, distributes these values, or portions of the values, among 
the p agents 3404-3407, receives intermediate results from the p agents 3408-3411, 
and returns computed values for the costate vector over a range of x values from 0 to 
Tk 3412 to the super. As discussed above, in currently preferred embodiments, the p 
module interpolates, rather than integrates, costate vectors, and thus does not need to 

20 distribute the problem in the manner illustrated in Figure 34. 

In the following discussions, for notational convenience, the u and r 
vectors are considered together as the u vector. In other words, in the subsequent 
discussion, the binding variable coefficients and barrier-function coefficients together 
compose a single control vector u. Recall that the iterative optimization strategy 

25 employed by embodiments of the present invention involves integrating — (r), 

expressed as follows: 

(~\- w (y( T )> u ( T )) 

dr K > p(y(r),u(T)) 

where 
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Hy(*)M*)H H (y(*)M*))TG(y(T),u{T)) 
p{y(T)>«(*))=( G (y( r )> u (T))) T 

5 In one embodiment, the integration of — (r) is carried out by a forward Euler 

dr 

integration method with variable step size At k : 

y<( Tk+ AT t ) = y i (T k ) + AT k ^-(T k ) 

Therefore, as is apparent in the above expressions, the y module, for each forward 
10 integration step, needs to invert the Hessian matrix in order to compute the w vector 
as the product of the inverse of the Hessian and the vector G Figure 35 graphically 
illustrates this computation. Inversion of the Hessian matrix is by far the most 
computationally intensive task involved in optimization. Therefore, the y module 
distributes inversion of the Hessian matrix among the y-module agents, relying on 
1 5 several assumptions for this distribution. 
Recall that: 

Figure 36 illustrates multiplication of the Hessian matrix by the w vector to produce 
20 the g vector in a distributed fashion. As shown in Figure 36, the G and w vectors, and 
the Hessian matrix, can be vertically partitioned into partitions i = 0 to 5. The first 
partition of the G vector 3602 is referred to as G°, the first vertical partition of the 
Hessian matrix 3604 is referred to as H°, and the first partition of the k> vector 3606 is 
referred to as w°. The Hessian matrix can additionally be partitioned horizontally into 
25 horizontal partitions j = 0 to 5. In Figure 36, the first horizontal partition of the 
Hessian matrix 3608 appears as a first vertical column within the Hessian matrix. 

Both vertical and horizontal partitioning of the Hessian matrix results 
in a block-like partitioning of the Hessian matrix into submatrices. For example, 
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Figure 36 shows the Hessian matrix partitioned into thirty-six submatrices, beginning 
with the left-most and upper-most submatrix J? 0,0 3610 and ending with the right- 
most and lower-most submatrix H 5 ' 5 3612. Note that, in the case of a diagonally 
dominant Hessian matrix, in which the diagonal block submatrices H 0,0 , H lt \ H 2 ' 2 , 
5 H 3 ' 3 , H 4 ' 4 , and H 5 5 have elements with magnitude much greater than the elements of 
other non-diagonal submatrices, then matrix multiplication of the Hessian matrix by 
the w vector can be approximated by matrix multiplication of w-vector partitions by 
the block-diagonal submatrix within the corresponding vertical partition of the 
Hessian matrix. For example, in a full matrix multiplication, the G-vector partition 
10 G° equals the sum of the matrix-multiplication products of the partitions of the h> 
vector with corresponding blocks of the first vertical partition of the Hessian matrix. 
Formerly, 

G' = H'w 
= Y 4 H ij w j 

7=0 

When the Hessian matrix is diagonally dominant, then the off-diagonal 
1 5 submatrices have relatively less magnitude and are of less importance in computation 
of a w-vector partition. This situation occurs when the Hessian matrix is partitioned 
in a way that removes, as much as possible, dependence between each partition and 
all other partitions. Therefore, a particular y-module agent can approximately 
calculate the current value of a ^-vector partition 3810 by inverting a corresponding 
20 Hessian block-diagonal submatrix 3808 and multiplying the large, bracketed term 
3812 in Figure 38 by the inverse of the block-diagonal Hessian submatrix 3808, using 
previously calculated values for the remaining H>-vector partitions 3814-3817. 

Figure 39 illustrates computation of a H>-vector partition w 1 by the 
second y-module agent. The second y-module agent 3902 obtains the current Hessian 
25 matrix block-diagonal submatrix H l ' } and current G-vector partition G 1 and most 
recently calculated values of the w-vector partitions w°, w 2 , w 3 , w 4 , and h> 5 from the y- 
module 2808, computes the inverse of the current Hessian block-diagonal submatrix 
H ]>] and then, using the computed inverse of the Hessian matrix, computes the current 
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approximate value of the w-vector partition w ] . The y-module agent 3902 can then, 
in a next computational step, furnish the approximate current value of the ^-vector 
partition to the y-module 2808 for distribution to the other y-module agents, and to 
receive their most recent computed values for the other n>-vector partitions h>°, w> 2 , h> j , 
5 w 4 , and w 5 in order to compute a next n>-vector partition value w 1 at the next x 
increment. Using this method, the y-module agents can synchronously compute the 
approximate value of n>-vector partitions at each step. 

Next, a control-flow diagram illustration of the super module is 
provided. Figure 40 is a control-flow diagram illustrating a highest conceptual level 

10 of the super module. In step 4002, the super module determines whether or not it is a 
replicate, or child thread, spawned or forked from a parent thread. If not, and the 
current thread is therefore the root thread, then in step 4004, the super module 
determines an initial point for the optimization and, in step 4006, initializes various 
computational entities, including the state vector, costate vector, and other sets of 

15 values, and also initializes a number of variables, including the Boolean variables 
"critical" and "precritical" to the Boolean value FALSE. Following initialization by a 
root thread, or as a result of flow of control from step 4002 for a replicate thread, the 
super module next chooses, in step 4008, a horizon T. Then, the super module 
invokes the function "integrate_w_and_y_to_T," in step 4010, which carries out 

20 successive inner-loop interactions of the optimization method, described above. 
Next, in step 4012, the super module invokes the routine "evaluate" in order to make 
decisions based on the inner-loop iterations carried out by the routine 
"integrate_w_and_Mo_T" in step 4010. The routine "evaluate" returns a Boolean 
value "continue" which the super module evaluates in step 4014, in order to 

25 determine whether or not to carry out yet another iteration of the outer loop of the 
optimization method, described above. If so, then the super module calls the routine 
"integrate p" 4016 to back calculate costate-vector values and then returns to step 
4008 to initiate another outer-loop iteration. If the super module decides not to 
continue outer-loop durations, then the super module invokes the routing "report" in 

30 step 40 1 8 and then returns. 
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Control-flow diagrams illustrating steps 4004, 4008, 4010, 4012, and 
4018 are next provided. The routine "integrate /?" called in step 4016, and the 
routines "integrate w" and "integrate y" called from the routine 
"integrate_w_and_jMo_T" are subsequently described mathematically, rather than 
5 using control-flow diagrams, in the interest of clarity and brevity. 

Figure 41 is a control-flow diagram for the routine "initial point" 
called in step 4004 of Figure 40. In order to determine an initial value for the state 
vector,^, the initial point routine, in step 4102, guesses at a value using any supplied 
hints or heuristics from users or from the mathematical modeling or preprocessing 
10 portions of the optimizer. Then, in the for-loop comprising steps 4104-4107, the 
initial point routine evaluates each constraint g, in order to determine whether or not 
the constraint, evaluated at the initially guessed state-vector pointy is greater than 0. 
If so, then the constraint is placed into partition / in step 4107. If not, then the 
constraint is placed in partition J in step 4106. Then, in step 4110, the initial point 
15 routine constructs an objective function F as follows: 

F = -5>, log( ft Crt)-5>> log(g y O0 + r y )-5>, log(2I) 

which is then minimized in step 4112 by a recursive call to the optimizer (Figure 30) 
to produce a modified initial state-vector value. The objective function Fcan be 
thought of as initially expanding the problem-domain boundary of the objective 

20 function F in order to envelope the initial value for the state vector, and the 
optimization of Fthen contracts the problem-domain boundary, forcing the initial 
point into the interior of the problem domain of the objective function F. In step 
4114, the initial point routine determines whether the modified state-vector value is 
feasible. A feasible state-vector value is an interior point of the problem domain with 

25 a reasonably high probability of being close enough to a local minimum for the 
optimization method to find a local minimum. If the modified state-vector value is 
feasible, then it is returned as the initial point for the optimization. Otherwise, 
control flows back to step 4102 for another pass at determining an initial state-vector 
value. In subsequent iterations, information gleaned from previous iterations may be 

30 employed in step 4102 in order to better guess an initial value for the state vector. 
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Figure 42 is a control-flow diagram for the routine "choose_next_T." 
In step 4201, the routine "choose_next_T" assigns local variable "T" to a minimum 
horizon constant. Then, in the for-loop, comprising steps 4203-4207, the routine 
"choose_next_T" calculates, for each y-module agent, a T value obtained by dividing 
5 a constant parameter by the minimum absolute eigenvalue of the block-diagonal 
Hessian submatrix associated with the y-module agent. If the computed T value for 
the y-module agent is greater than the contents of the variable "T," as detected in step 
4205, then the variable "T" is assigned to the T value calculated for the y-module 
agent. Next, in step 4208, the routine "choose_next_T" initializes local variables "n" 

10 to a maximum number of increments, initializes local variable "a" to an increment 
coefficient, and initializes a loop-control variable "j" to 0. Then, in steps 4210-4212, 
the routine "choose_next_T" increases the value of the variable "T" in step 4212 until 
the derivative of the costate vector with respect to T evaluated for the current contents 
of the variable "T" approaches 0, as detected in step 4210. The loop may also 

1 5 terminate when the number of iterations is equal to the contents of the local variable 

"n." When ~~ = 0, then computation of the costate vectors over the continuous 

iteration range 0 to T can be more easily calculated. In alternative embodiments of 
the present invention, the costate vectors can be computed using a polynomial 
approximation, rather than employing complex integration. 

20 Figure 43 is a control-flow diagram of the routine "evaluate," called in 

step 4012 of Figure 40. In step 4302, the routine of "evaluate" sets the Boolean flags 
"atjninimum" and "continue" to the Boolean value FALSE. Next, in step 4304, the 
routine evaluate determines whether the computational resources allocated for the 
optimization have been exhausted. This may include continuous computation for 

25 longer than a maximum specified time, the usage of processing resources, memory, or 
a combination of computational resources greater than a specified maximum amount 
of resources to be used, the occurrence of error conditions, or other types of 
predetermined termination criteria. If the computational resources are exhausted, 
then the routine evaluate returns. Otherwise, the routine evaluate, in step 4306, 

30 determines whether the optimization thread is currently critical. An optimization 
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thread is critical if all >>-module agents are currently in Morse mode, and have been so 
for some predetermined number of iterations. Morse mode is described in detail, 
below. If the current thread is not critical, then the routine "evaluate" sets the 
Boolean flag "continue" to TRUE, in step 4308, and then returns. Otherwise, in step 
5 4310, the routine "evaluate" determines whether the thread, currently determined to 
be critical, has reached a local minimum. If so, then the Boolean flag "at_minimum" 
is set to TRUE, in step 4312, and the routine "evaluate" returns. If the optimization 
thread is not at a local minimum, but is critical, then the routine "evaluate," in step 
4314, determines a number of replicant threads to spawn, and also initial directions or 

10 starting points for the replicant threads. Spawning of replicant threads is illustrated in 
Figure 27. In the for-loop of steps 4316-4319, the determined number of replicant 
threads are spawned, each replicant thread linked has a child to the current 
optimization thread. Then, in step 4320, the routine "evaluate" may attempt to 
perturb the current optimization trajectory in order to move the optimization 

15 trajectory away from the critical point and continue descending towards a local 
minimum. In step 4322, the routine "evaluate" sets the Boolean flag "continue" to 
the Boolean value TRUE that then returns. 

Figure 44 is a control-flow diagram for the routine "report." This 
routine is invoked in step 4018 of Figure 40. In step 4402, the routine "report" 

20 determines whether the current optimization thread has a parent. If not, and the 
current optimization thread is therefore the root optimization thread, then the routine 
"report" determines, in step 4404, whether there are any children threads still running 
that are linked to the current, root optimization thread. If so, then the current, root 
optimization thread waits, in step 4406, for the children optimization threads to finish 

25 execution. Next, in step 4408, the routine "report" sets local variable "min" to the 
current value of the objective function computed with the current state-vector value, 
sets the local variable "location" to the current state-vector value, and sets the local 
variable "at min" to the Boolean value flag "atminimum" returned by the routine 
"evaluate" for the current optimization thread. In step 4410, the routine "report" 

30 determines whether any children threads have terminated and reported back to the 
current optimization thread. If so, then in the for-loop of steps 4412-4415, the current 
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optimization thread evaluates the report generated by each terminated child in order 
to determine the lowest value for the objective function reported by any child and 
achieved by the currently executing optimization thread. Next, in step 4416, the 
routine "report" determines whether there are any children still running. If so, then in 
5 the /or-loop of steps 4418-4420, the routine "report" links each child to the parent of 
the currently executing thread. Note that the test in step 4416 fails in the case that the 
currently executing optimization thread is the root thread. If the currently executing 
optimization thread is the root thread, as determined in step 4422, then the routine 
"report" returns the minimum objective function value, state vector, and indication of 

10 whether the optimization is currently at a minimum in step 4424. Otherwise, the 
routine "report" reports these same values back to its parent thread, in step 4426. 

Figure 45 is a control-flow diagram for the routine 
"integrate_w_and_y_to_T," called in step 4010 of Figure 40. In the for-loop of steps 
4502-4507, the routine "integrate_w_and_y_to_T" successively calls the routines 

15 "integrate w" and "integrate y" in steps 4503 and 4504, in order to complete one 
inner-loop execution for the optimization method. The successive iterations continue 
until either the optimization thread goes critical, as detected in step 4505, until the 
value Tk is greater than or equal to T, or until some other stopping criteria are 
detected, such as exhaustion of computational resources or other such stopping 

20 criteria in step 4507. 

Figure 46 is a control-flow diagram for the routine "integrate^." This 
routine describes the operation of the y module (2808 in Figure 28). In step 4602, the 
routine "integrate /' obtains symbolic representations of the Hessian matrix, the 
vector G, and other data needed by the routine "integrate /' from the super module. 

25 Next, in step 4604, the routine "integrate /' computes numerical values for the 
Hessian matrix and vector G*. In step 4606, the routine "integrate yp distributes 
partitions of the computed Hessian and G vector to each ^-module agent. The y- 
module agents compute their associated »v-vector partition values and distributed row 
values, which are collected by the routine "integrate y in step 4608. In step 4610, the 

30 routine "integrate y" determines whether any agents are currently in Morse mode. 
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The term "Morse mode" refers to an agent that currently needs to 
compute its w-vector partition based on the Morse quadratic approximation of the 
original objective function associated with the y agent method since the agent is 
currently near a critical point with respect to the agent's partition of the vector G. If 
5 any agents are in Morse mode, then, in step 4612, the routine "integrate /' directs the 
Morse-mode agents i to recompute p l and W and collects the recomputed />' and W 
from each Morse-mode agent /. Next, in step 4614, the routine "integrate /' 
determines whether all agents are currently in Morse mode. If so, then, in step 4616, 
the routine "integrate /' determines whether the local variable "precritical" currently 

10 has the Boolean value TRUE. If so, then the agents were previously in Morse mode 
in the previously executed iteration, and so the local variable "criticalcount" is 
incremented, in step 4618. In step 4619, the routine "integrate /' compares the 
contents of the local variable "critical_count" to a threshold constant and, when the 
contents of the local variable "critical count" equals or exceeds the critical constant, 

15 the routine "integrate y* concludes that the current optimization thread has remained 
critical for a sufficient number of iterations to declare the optimization thread to be 
critical by setting the Boolean flag "critical" to Boolean value TRUE in step 4620. 
However, if the contents of the variable "critical_count" do not justify declaring the 
current optimization thread to be critical, then control flows to step 4624. If, back in 

20 step 4616, the routine "integrate/' determines that the variable "precritical" did not 
have the Boolean value TRUE, then, in step 4622, the variable "precritical" is set to 
true and the variable "critical_count" is initialized to 0 before control flows to step 
4624. In step 4624, the routine "integrate y computes global p and iv-vector values, 

and it computes from them the value of — . 

dr 

25 Next, in step 4625, the routine "integrate /' determines a value Ax for 

the current Euler forward integration step. If the computation of Ax succeeds, as 
determined in step 4626, then the routine "integrate y carries out the forward 
integration step, in step 4628, updates the binding variable coefficients, if necessary, 
in step 4630, and returns. If the computation of Ax fails, as detected in step 4626, 

30 then the routine "integrate y returns with an error. 
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The forward Euler integration technique has been discusses above, but 
is presented formally below. The global value for p is assembled from the partitioned 
p values computed by y-module agents in a current inner-loop step as follows: 

n—\ 

pW r *)»»( r *))=Zp , W r *)» i< ( r *)) 

1=0 

5 Next, the values for the partitions of the w vector are collected from the ^-module 
agents and assembled to produce a w vector. Then, a continuous iteration increment 
Ax is computed, by a method to be described below. Finally, the state vector is 
propagated from a current continuous iteration value r* to a next continuous iteration 
value Tk+dr as follows: 

10 y'(T k+ AT k ) = y i (T k ) + AT k ^-(r k ) 

Finally, the continuous iteration value r* is updated as follows: 

n =r A +Ar 4 

Note that a jy-module agent is considered to be critical with respect to the ^-module's 
vector-G partition when the absolute value of the^-module agent's local p, p\ falls 

1 5 below some constant parameter. The computation of a w-vector partition W and the 
y-module-agent-local p, p\ is described below. 

Figure 47 is a control-flow diagram for the routine "compute Ax," 
called in step 4625 of Figure 46. In step 4702, a guess, or estimate, Ax', is made for 
Ax. This guess or estimate may be obtained from a constant parameter, or may be, in 

20 more elaborate embodiments, derived by empirical methods or heuristics from prior 
computed values of Ax and from recent history of state-vector values, control-vector 
values, and other such information. In the described embodiment, for forward Euler 
integration steps, when the control vector has an approximate constant value, then F 
decreases by one at each step: 

25 
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dF 

When some of the ^-module agents are in Morse mode, may depart from the 

dr 

dF 

value However, if the value of ever rises above 0, as detected in step 4704 

dr 

of the routine "compute Ax," then an error condition obtains since — (r k ) is not a 

dr 

5 descent direction and therefore the forward Euler integration method cannot produce 
a feasible trajectory y(x). Otherwise, the estimate of At, At 1 , is inserted into the 
forward Euler integration step equation 

y(T k +Ar') = y(T k ) + AT'^(T k ) 



10 and the state vector is propagated to y(r k + At'). This computed next state-vector value 
is checked, in step 4706, for feasibility. If the state-vector value is not feasible, then 
the estimate of Ax, Ax', is decreased, in step 4708, and a new forward state-vector 
value is again computed and checked for feasibility in step 4706. Once a feasible 
forward state-vector is found with a corresponding Ax estimate, Ax', then, in step 

15 4710, the Armigo condition 

F(^(r i+ Ar') )M (r t ))- J F( > ;( rt ) jM ( ri ))< 0 .A T 'G r ( r (r i ) ) «(r,))^(r t ) 



is evaluated with respect to the estimate of Ax, Ax 1 . If the Armigo condition is 
satisfied, then the estimate of Ax, Ax', is returned as the new computed Ax. 
20 Otherwise, Ax is approximated. Several different approximation methods are 
available. In the case that a previous feasible estimate of Ax, Ax A , is available, then a 
quadratic approximation to the objective function F around x k can be made, as 
follows: 

F(y(r k + Ar),w(rJ)«F(Ar) = aAr 2 +Mr + c 

25 

Then, an estimate for Ax is obtained as follows: 
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F(0) = F(y(r k ),u(v k )) = c 
P(At a ) = F[y (r k + Ar A ),« (r t )) = a(Az A f + b (At a ) + c 



(0) = {G(y(T k ),u(T k ))f%(T k ) = b 



dF 



= F(0) 



6 = (G(>,(r,), M (r,))) r ^(r t ) 



F(A^)-F(0)-A^(G( 7 (r,),«(r t )))^(r t ) 

# = £±± 



A At along the descent direction that minimizes F can be approximated by finding a 
Ax that minimizes F , which occurs when: 

15 

dF (At) 
,. v / =2aAr + b = 0 
d(Av) 

Therefore: 

Ar' = -— 
2a 



20 When two previous feasible estimates for Ax for the current x k are known, then a 
cubic approximation method may be employed. However Ax is approximated in step 
4712, the new approximate Ax is used to recompute a forward value for the state 
vector, which is then tested for feasibility in step 4706. 
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As discussed above, when a ^-module agent is near a critical point and 
is therefore in Morse mode, the standard forward Euler integration procedure 
discussed above is not useful in finding an optimal optimization trajectory y*. In 
general in the neighborhood of a critical point, the Hessian block-diagonal submatrix 
5 Z/ 11 corresponding to ajy-module agent / may have negative eigenvalues, and therefore 
the computed direction might not be a descent direction. However, the Hessian 
submatrix H* 1 can be modified in such a way as to produce a positive definite Hessian 
submatrix H' 1 in order to carry out a Morse-mode forward integration of the state 
vector. 

10 The derivation of the Morse-mode approach to forward integration of 

the state vector is based on the following three assumptions: (1) the objective 
function F(>>(r),w(r)) is roughly constant with respect to u(x)\ (2) the Hessian 

matrix (^(r),w(r)) is diagonally dominant, implying that thejy-module agents are 

partitioned for relative independence; and (3) the state vector y l (x) changes very little 
15 while ^-module agent / is in Morse mode. Assume that ^-module agent i enters 
Morse mode, or local criticality, at r k = f ' , as determined by the value of 
/?'0(r),w(r)) calculated by the standard method. Then according to the first and 
third assumptions, above, the objective function F may be approximated about the 
domain point (y 0 (T k )... 9 y°(T i ) 9 ... 9 y n ~\T k ),u(? i )) by a Taylor approximation as 
20 follows: 
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F(Hr),«(r)) = F(/(r t )...,y(f'),.^-'(r 4 ),«(f')) 
+(G'(/(T t ).,y(f)...o^(r 4 ).«(f))) r (y(r)-y(f)) 

+|(G'(/(rj... ) y(f'),..y- l (r t ) > »(f'))) r (yw-y(-')) 

4(^( r )-^( f O) r ^1^(ro...,y(fO,..r- i (rj, M (fO)(y(^)-y(^0) 
+(y(o-y(#')) r z« i '(y , (^)-y(# , ).-y rt ('t).«(#'))(yw-y(r.)) 

/=0 

/*/ 

+{S(y«V(rO) r i^(y(rJ...y(f') I ...y^(r t ),«(#'))(y , (r)-y 
^(|y(r)-^(f)|f) 
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Dropping the third-order term, the objective function can be approximated as: 



F(^(r) ) «(r)) = F(/(r t )... s y(f'),..y"- 1 (n),«(f')) 

+(G'(y , (r 4 )..,y(f),...y rt (r t ) i «(f))) r (y(r)-y(f)) 
+|(G'(/(rj...,y(fO,..y-'(r 1 ),«(^))) r (y(^)-y(^)) 

+(y(O-y(^0) r Z^1^(rj...,y(f0,-..y-'(r t ), M (f'))(y(r)-y(r i )) 

/=0 

/*» 

+ i$y (0 V K)) r ^ (y (r,)...,y (f ),..y-(r t ), M (f'))(y' (r)-y' (rj) 

The block-diagonal Hessian submatrix ft 1 around the domain 
5 point (y 0 (T k )...,y 0 (? '),..., y n ~ ] (r k ),u(f')) can then be transformed to a positive definite 
matrix using a weak Morse transform, as follows. 

First, as shown in Figure 48, the Hessian matrix can be decomposed by 
spectral decomposition as follows: 

HP = PA 

10 

HPa^KPa *ra=0 N-l 



U ••• o 1 



A = 



'■N-l J 



15 



H = PAP r 
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The NxN matrix P 4802 comprises eigenvectors for the Hessian matrix as columns, 
and the NxN diagonal matrix A 4804 includes the corresponding Eigenvalues of the 
Hessian matrix. 

5 One approach to transforming the Hessian is to create a new matrix 

P = [p 0 , p v p n _ x ] in place of the matrix P, where the column vectors are defined as: 



p A Pa 

a lo 



if K >o] 

if/l„ <0 



>for a = 0, N - 1 



10 



In addition, a new matrix A = 



0 



0 



is defined with diagonal entries: 



A "=o vcJ*"- 0 N -' 



Using the previous results, a modified Hessian matrix that is positive definite can be 
obtained from the modified P matrix and the modified A matrices as follows: 
15 H = PAP T 



Furthermore, a reverse for the modified Hessian matrix 
H can be computed as follows: 

H R =PA«P T 

20 A reverse is needed because the transformed Hessian may be singular. 

A second approach to transforming the Hessian is known as the 
Strong-Morse Transform, and is used in the Morse-mode forward integration step. In 
this method, a new modified A matrix is produced as follows: 
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4=K|for a=0,...,N-l 

Then, a modified Hessian is computed as follows: 

// = / > AP r 

5 

This transformed Hessian is non-singular and thus the inverse of the transformed 
Hessian is well defined, as follows: 

ir x = pa.- 1 p t 
» -i i 

10 A =4 

aa j 

A, a 

In alternative embodiments, the weak Morse transform can be used for points that are 
outside of the Morse mode, points that are far from a critical point. 

Using the transformed Hessian, H u \ the approximation of the objective function 
15 becomes: 
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+(G , (/(r,)..,y(f)....y^(r t ).«(f))) r (y(r)-y(r')) 

+|(G'(/(n)--.y(-0'-y- , (-J."(^0)r(^( 7 )-y(-0) 
+(yw-y(^)) r g^(y , (^)---.y(#')...-y M (r*).«(#'))(y(r)-y(T 4 )) 



Taking the derivative of this approximate objective function with respect to the i 
partition of the state vector y* (r) yields the ith partition of the vector G, as follows: 



= G'(y\r k )...,y'(T%.y-\T k ),u(r')) 

+^(y , (r t )...,y(f)....y-'(r 1 ),«(f))(y(r)-y(f)) 

+z^1/(^)--.y(^)»-y-'(^)."(^))(^w-y(^)) 



5 The <* partition of the vector G evaluated at x k is then: 
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0-(,(r),.(r))-l ^ 1 » \ 

T ~ T k 

= G'(/(r i )...,y(f'),...^- , (r t ) > «(f')) 

+ ^(/(rj...,y(fO,.^- , (rO, M (f))(y(r t )-y(f')) 

The Hessian submatrices are therefore given by: 

mi < w « I WM-An-jTM-WP'-J 

[^(/(r t )...,y(f'),..y- , (r), M (f'))otherwise 



5 Again, the off-diagonal Hessian submatrices are not modified by the Strong-Morse 
Transform since the Hessian is diagonally dominant, and the off-diagonal submatrices 
contribute very little to the approximation of the objective function. Using 
assumption 2, above, the block-diagonal Hessian submatrix at z k is approximately 
equal to the block-diagonal Hessian submatrix at f ' 

10 H»(y°(T k )...,/(?),..y-*(T k ^^ 

Similarly the third assumption, above, implies that 
G'(/(rJ...,y(fO,^ 

15 

Therefore, the vector G of the approximated objective function 
F(y(?)>u(r)) evaluated at r k is well approximated as follows: 

20 ^WT*).«(r*))-^Wr0.«(r0) + ^(y(# , ).«(^))(y(r0-y(f')) 
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and the Hessian submatrix for the approximated objective function evaluated at \ k can 
be approximated as 

- / / ^ / ,v f H"{y{z i ),u{r'))\ii = j 
\H J {y{r k ),u(T k ) J otherwise 

These approximations then allow for approximate calculation of the ^-vector 
5 partition and the local p for agent i, w l and p\ as follows: 



^W^).«K))-Z^(j'(«*).«(^))^(j'(«-*-.).«(«- m )) 



y=0 



p'(^(r t ),«(r t )) = (G'(y(r t ),«(r t ))) r V(j;(r 4 ),«(r t )) 

i=0 

10 

In step 4630 of Figure 46, the binding variable coefficients may be 
adjusted. Although the binding variable coefficients might be continuously adjusted, 
as implied by Figure 19, in at least one embodiment of the present invention, the 
adjustment is binary and occurs when 
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y N - C{y)+ J] r a s a <u M + ipm__R_switch_of f set 
V «=o ) 



where y N is the envelope variable z. 

At the point when the binding variable coefficients are adjusted, they are adjusted to 
some larger constant value from their original smaller, constant value. At this point, 
the state vector needs also to be adjusted, so that adjustment of the binding variable 
coefficient does not force the state vector to an infeasible point. The state vector is 
adjusted as follows: 



where the parameter "ipm_c_offset2" has a value close to one. 

The w_and_r module needs to compute a next value for the control 
vector u during each iteration of the inner loop. This computation is very similar to 
computation of y vector in jy-module agents. As discussed above, the iterative 
equation for calculation of a next control vector u is: 



approximation of the Hamiltonian for each iteration. Since we can find an explicit 
expression for the minimum we do not necessarily need the parameter a. The u 
computation only takes one step. 



approximation around u k , given a fixed yk and using a second order Taylor 
expansion can be computed as follows: 




For the formulation of the static problem we use a convex quadratic 



Let u k = 



w| r ^be a point of the function li,(y,u,p). A second order 
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H\ «H\ +T| ( u - Uk ) + L(u-u k ) T — 



yk>»k>Pk 



d 

where is defined as . This is a quadratic function, but it is not necessarily 

du 

dW 

convex since is not necessarily positive definite. A Morse-transform can be used 

du 

dW d^ 

to transform to a positive definite convex approximation , as was done above 

du du 

5 for the Hessian matrix. 

The Hamiltonian matrix can then be approximated as: 



Hi *U\ +W\ (u-u k ) + -tu-u k ) 7 

l^.w.p* *yk> u k,Pk hk,"k,Pk y k/ 9 V 7 



du 



yk >"k >Pk 



Defining Au = u-u k , 



■^.Aii.ft hk^kiPk *yk*k>Pk 2 du 



Au 



yk>"k*Pk 



10 In the following discussion, the dependence on w*, and /?* is not explicitly shown. 
The minimum of the above function can be found by setting the gradient of the 
objective function equal to 0 and solving for Aw. Since 



dH T d"¥ 

dAu du 



y*k,p 



then 



15 



where the superscript R denotes the reverse, described above. Defining: 
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then 



Au = -a 



5 This allows for iterative update of the control: 



=u k +Au 



However, solving for *F and is a computationally involved process. The 

du 



Hamiltonian is expressed as follows: 



10 



= 1_£z2 



Using the chain rule and the fact that p is treated as a constant during the control 
1 5 computation: 



« f i i 5 P 



for a = 0, . . M - 1 and 



cW a _ d 2 H 



du p du a du p 
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N-i 



1 d w c 1 dw c dp 1 dw c dp I d 2 p 2 dp dp 



w ' 1 w 

P 9u a du fi P 2 d" a du p P 2 d" a P 2 ( d» a d"p P 3 ( du a du 



for <x,P = 0, . . ., M - 1 
By definition: 



5 for a = 0, . . N-l, the first derivative of which, with respect to w, is: 



f dH~ l 



du 



du. c a < du 



f* J 



fora = 0, AM and/? = 0, ...,M-1. 

The inverse of the Hessian matrix H~ 1 is defined as: 
10 H~ X H = I 

and can be written in terms of its elements as: 



<r=o 



' H 'al H Cn = S arj 



1 5 for a,p = 0, . . N - 1 and where: 



arj (0 if otherwise 



Taking the derivative with respect to u on both sides, using the product rule: 



20 



N-l 

<r=o 



du ( " " 



du. 



= 0 



r J 
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for ot,T| = 0, . . N - 1 and y =0, . . .,M-1 . Rearranging terms: 

£=! dH 'ac S=J 1 dH c n 

<-=o ou r n=0 ou y 

5 

and multiplying both sides by H~ ] : 



N-\ /V-l r)fj~ l N-l N-\ p)JJ 



10 the expression can be simplified to: 



"=! dH~ l , dH~\ 

y — = — ^> 

^ 5w r hP du y 



providing an explicit expression for the elements of H 1 



15 - = ~YY H aC -Kl 



Substituting the expression for H 1 into 

N-\ 

=z 



3w_ %=!( dH:\ _ dG, ^ 



K du p < «du tJ 



20 provides: 
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dw. 



du 



0 <T=0 



dHZl 



du p ( a€ du a 



P J 

N-\( ( N-l N-l F)JJ ~\ 



^ 77=0 co=0 uu p j 



p J 



N-l N-l N-l p)H N-l fifZ 



N-l N-l 



, =0 du fi 
SO, 



77=0 eo=0 n dllp 

r 



rj=0 



77 o 



fora = 0, ...,N-\ andy!? = 0, ...,M-1. 
The second derivative is 



du p du y 7=0 



5G, 



arj 



y t=i du p du r j 



Note that the expression is simplified by the fact that the second derivative of the 
gradient of F and Hessian matrix with respect to u is 0 since these terms are linear in 
u. Substituting 

the expression for H~ l into the above expression provides: 



du R du„ r=a I 0=o du y 

7=0 

du 



N-l ( N-l fi]-[ rs \ 

?j=0 ^=0 "Up ou r j 
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This can be simplified as follows: 



9«^aw y f^o [££ du y du p 



-Yh-1 Y- 

v=o ^S3> dup du y 



N-l 



\ 



c=o ^=o ou r j c=0 ^ aw^ aw y ; 



£-1 ( N-l Qff 



CO vw e 



= -Yh~ 1 c Y 

c=0 [po du r dup dup du Y j 



Next the denominator, p, is defined as 



The first and second derivatives with respect to u are, using the product rule, 



dp 



du 



N-l 

=z 

a (=0 



3G, 



<- w +G ^£ 
5m„ c c du 



a J 



fora = 0, ...,M-1 



and 



d 2 P 



du a du p (=0 



dG, dw, dG, dw, d 2 w r ^ 

— ( - — ^ + — ( - — - + G C £- 

K du a du p du p du a du a dup } 



for a, p = 0, ...,M-1 
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Next, a formal distributed algorithm is discussed. Multiplying both sides of: 



by: 



du 



results in: 



du ' 



Expressing the relationship in terms of the r* agent 



10 for i = 0, . . ., M- 1 . Rearranging terms: 



Using the notation: 



du- du s * 



s=0 



du r 
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Multiplying both sides by the inverse of provides: 



d*¥ r 1 - 

du r 



results in: 



5 Thus, the control can be updated by: 



Au r = -<f> r 



The definition for the Hamiltonian is now written in terms of agents: 



n-l N'-\ 



_ i ;=0 f=0 



10 Then 



\pr _££L__v v 



du 



a i=0 <=0 



1 dw> e l_ ,. dp 

P du r a P 2W "du r a 



and 



A\JJ r n-l N'-l 

ou p /=o c=o 

1 d 2 w^ l m£ dp 1 dw£ dp 1 , d 2 p 



dp dp 



\\; 1 — i _ 

_p du r a du; p 2 du; du r a p 2 du 5 a du; P 2 * du r a du; P 3 * du r a du 



p J 



Considering that: 
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for a = 0,...,n l and that: 



n-l N J - 



5 each element a of the w vector for agent I can be expressed as: 



aw 



/» <T=o 



v ^ G' +(//") ^ 

du r a c v '<* a«: 



Substituting the definition of H 1 in the previous expression: 



10 



N' -I N' -1 



^ =0 II ^=0 o=0 

Rearranging terms: 



PhJ N*-\ , fiCi* H'-ltf-l , nr -, 

— = Y ("") — - X H ("T — ^ 



15 



du r p du r p 



Since the second derivatives of G and H with respect to u are 0: 



du r fi du r f=0 



du 



However: 



f§L-**£!Sv Wr fy 

k a«; £s a W/ , "J ^ a*; a*; 



^ a 2 G^ 
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aV 



N'- 



r 



N'-l N'-l 



rHV x dH'^dw^ 



— = Y -Y.Y ( H ")~ l 



k=0 0=0 



* ^ awi a«! 



o>-0 



f d 2 G, ^ 

y du r p du s rj 



Separating terms and changing indices of the last two terms: 



d 2 < 

du;du; 



V 

_y y (h"\ 0 



-Y y (h"Y «° 9 



JV'-l . w'-l fill" p^J \ N'-l 

<o=0 S fi>=0 OU p OU y J £=Q * 



r d 2 & ( > 



AT=0 



a//" dw' 



N'-i , A/'-l 



k6 



k=0 



e=0 y du; du; du; du ! r j 



9=0 du; du; I 

AT=0 



K du;du r) 



/V -1 



d 2 G' K 

y du r /u; 



Therefore: 



r)Ci' rid' " -l A,7 - ) 



15 



— + — L 

v a«; < «du' PJ 



and using the fact that the second derivatives of G and H with respect to u are 0: 



paf^i n-l N j -\f 

Finally: 



,=o <;=\-o 



dH 1 ^ dw J , dH'l dwi 



a«; aw; 



aV 



and its derivatives are: 
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du r a du; 




«^sa«; a^ + a^a< 




The propagation of the costate variables is different from that of state 



and control variables, which are tightly coupled and proceed together forward in time. 
In contrast, the costate variables are propagated backwards in time, independently of 
the integration of the state and control variables. In order to do this, a database of 
values calculated in the y module and w_and_r module is needed. Since the 
10 differential equation of p is linear in p, its solution can be approximated by an 
interpolation method rather by integration. For many applications the value of the 
vector p is almost constant with the exception of the transient around the convergence 
horizon. Also the matrix ,4 in the equation for p does not depend on p. 



15 calculated by the y module and w_and_r module is created for use by the p module. 
The terminal conditions p(T) are computed first, and then earlier values of p are 
calculated in a loop which interpolates backwards to the beginning of the y 
computation. 



T k-\ ~ T k where t^i < where the matrix A is assumed to be nearly constant. Then, 
integrating backwards the solution at x^i given t* is 



To compute p 9 a set of the node values of y and u trajectories 




20 




For practicality, we divide the problem up into increments of size 



25 



P ( r *-i ) = p( T k) ex p{A {r K _ x - r K )) 
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The exponential of this equation is solved by the approximation 

exp(A(r /c _ 1 -t k )) = X V Vt 

where N is the degree of the approximation. 

5 Repair 

So far, the discussion has focused on specific computational methods 
for generic optimization. However, an additional aspect of real-world optimization 
concerns repeated optimization over real periods of time. For example, an enterprise, 

10 such as a manufacturing facility, may employ optimization techniques for optimizing 
all the various aspects of running the manufacturing plant, including scheduling of 
employees, ordering of parts and raw materials, distribution of parts, raw materials, 
and intermediate manufactured articles within the manufacturing plant, and delivery 
of product to customers, among other things. In general, in the extremely high 

1 5 dimensional problem domains for such enterprise optimizations, and because of the 
many uncontrolled events and changes that arise during operation of a manufacturing 
plant, optimization may need to be carried out repeatedly, over time, in order to 
readjust the control of the manufacturing plant to respond to changing conditions. 
Similar considerations apply to many other computational tasks, in addition to 

20 optimization. For example, modeling of a complex system may require intermediate 
computation of system state in order to respond to various randomly-generated events 
and perturbations to the system, introduced so that the model reflects real time 
response of the system to a changing environment. 

Figure 49 illustrates this repeated optimization strategy. In Figure 49, 

25 and in subsequent Figures 50 and 52-56, the horizontal axis 4902 corresponds to time. 
In Figure 49, a near-optimal state function for a system, including state variables and 
control variables, is computed via optimization and/or modeling techniques, at each 
point in time represented in Figure 49 by a vertical dashed line, such as vertical 
dashed line 4904. The computed state functions are represented in Figure 49 as 

30 curved segments, such as curved segment 4906, that begins at a point in time 4908 at 



76 



which computation of the state function begins and extends to another point in time 
4910. Of course, computed state functions are, in real-world applications, complex 
control trajectories in hyper-dimensional spaces, but two-dimensional control 
trajectories, such as curved segment 4906, are used for clarity of illustration. The 
5 diagonal, dashed lines, such as diagonal dash line 4912, represents the computation of 
the state function via integration with respect to a continuous iteration variable x, 
described above, and takes some amount of time 4914 related to the computation 
horizon T for the optimization technique by which the state function is calculated. 
When a new state function is calculated, as, for example, new state function 4916 

10 applicable to the range of time beginning at time 4918 and ending at time 4920, the 
new state function likely relatively closely overlaps, but does not coincide with the 
previous state function, in the present case state function 4906. Thus, there is a 
discrepancy 4922 or discontinuity between successively calculated state functions. 
This discontinuity arises from a change in conditions and environment of the system 

15 between the time the previous state function was calculated and the time that the 
current state function calculation began, in this case the events and changes that occur 
in time interval beginning at time 4908 and ending at time 4918. 

Unfortunately, these discontinuities represent disturbances in the 
system. For example, in the manufacturing plant, a discontinuity may reflect a 

20 change in scheduling of workers, a change in departure times and destinations for 
delivery trucks, and other such changes. However, complex systems tend to have 
inertia, and major disturbances are difficult to overcome, may take finite amounts of 
time and resources to overcome that themselves change the environment of the 
system to a degree that the previously calculated optimal state function is no longer 

25 optimal, and may cause ripples and cascading effects that further perturb a complex 
system. 

One possible approach to ameliorating the discontinuities arising from 
repeated optimization is to shorten the time between optimizations. Figure 50 
illustrates a shortened interval between optimizations for the state function 
30 calculations illustrated in Figure 49. By more efficiently computing the state 
functions, represented in Figure 50 by increased slope of the dashed lines, such as 
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dashed line 5002, representing state function calculation, the optimizations may be 
carried out in a shorter amount of time AT 5004. However, the effect of doing so 
increases the frequency of disturbance, and the increased frequency of disturbance 
may itself create rippling and cascading effects that eventually lead to greater and 
5 greater instabilities. In certain cases, even greater discrepancies, such as discrepancy 
5006, may arise as a result of underdamping of the effects of the occurrence of certain 
kind of events in very short time intervals. 

Figure 51 illustrates an analogy for a second consequence arising from 
repeated optimizations. Note that, in Figures 49-50, it is assumed that the state 

10 functions are calculated anew without regard for the previously calculated state 
functions, other than at the point where the previously calculated state function is 
abandoned and a new state function begins to control the system. However, failure to 
take into account the history of previous state function calculations may result in state 
function solutions that are exceedingly non-optimal over longer periods of time. As 

1 5 an analogy, consider the task of traversing terrain represented by a simple topographic 
map, in Figure 51, from point A 5102 to point B 5104. A first path 5106 from point 
A towards point B represents a path that might be followed by an individual relying 
only on a compass, who frequently readjusts his or her trajectory based only on an 
instantaneous compass reading without even waiting for the compass needle motion 

20 to settle. Such a path represents an exceedingly non-optimal solution, in time and 
distance, for traversing the terrain between points A and B. 

A second path 5108 is a path that might be taken by an individual who 
recognizes that the compass must be allowed to settle, before reading, who recognizes 
that a compass may be perturbed by local conditions, such as an outcropping of 

25 magnetite in a cliff face, and who recognizes that if an obstacle is encountered, it may 
be prudent to follow that obstacle rather than relying on the compass for direction. 
Obviously, the second path 5108 represents a far more efficient and effective strategy 
for traversing the terrain from point A to point B. 

A third path 5110 represents a path taken by an individual capable of 

30 reading a topographic math, properly using a compass, who recognizes that a bridge 
across a river is a worthwhile intermediate point to aim for, and who recognizes that it 
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is far easier to walk along a path through the forest than to attempt to blaze a new 
trail. In the analogy of Figure 51, the individual who can rely on past information, 
rely on information in addition to a simple, limited current state, and who can project 
future solutions based on additional information is the individual who is most likely 
5 to find an optimal path between point A and point B. Returning to the repeated 
optimization problem illustrated in Figures 49-50, it is desirable to find a way to 
employ previously calculated state functions in order to employ more information 
over longer periods of time, rather than to simply consider local, instantaneous state 
information. 

10 As a result of the considerations described above with reference to 

Figures 49-51, an approach to a more or less continuous optimization, over time, 
using historical information and updated projections, represents an additional feature 
of embodiments of the present invention directed towards optimization, referred to as 
a "repair-based continuous optimization approach" ("RCO"). Figure 52 illustrates the 

15 concept of a sliding window. As shown in Figure 52, a prior window 5202 beginning 
at time tl-AT 5204 and ending at time tl+T-AT 5206, of time-length T 9 is used as a 
frame of reference for computing a state function. Since that time, an amount of time 
AT 5208 has elapsed. Therefore, a new, current window 5210 of length T in time 
provides the frame of reference at the current time // 5212. 

20 As shown in Figure 53, the state function calculated in the prior 

window 5302 extends from time ti - AT 5204 to time ti - AT+T 5206. Figure 54 
illustrates extension of the previously calculated state function to extend over the 
current window. As shown in Figure 54, the previously calculated state function 
5303 is extended 5402 to the end of the current window at time t } + T 5404. Figure 

25 55 illustrates consideration, at the current point of time, of past events and 
environmental changes that may result in the previously calculated state function 
being out-of-date, and not optimal, at the current point of time, as well as an approach 
to incorporating information from the previous state function in finding a new, 
optimal state function for the current frame of reference represented by the current 

30 window. 
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In Figure 55, vertical lines 5502-5505 represent events that occurred 
during the time of calculation of the state function, ti - AT, and the current time, t\. 
Such events may alter the state and environment of the system in such a way that the 
previously calculated state function is no longer optimal. For example, in the 
5 manufacturing plant, a large number of employees may have decided to go on strike, 
or a flu virus may have greatly increased the number of workers unavailable for being 
scheduled to perform tasks within the manufacturing plant. Therefore, previously 
calculated optimal scheduling may be exceedingly suboptimal in the face of a 
shortage of workers. As shown in Figure 55, various projections of time-dependent 

10 parameters on which the state function was calculated may have changed at time t\. 
The change in these projections is represented by the dashed line 5506. For example, 
in the context of the manufacturing plant, a major customer may have indicated an 
immediate need for twice as many products as the customer normally buys, or the 
arrival of a new product produced by a competitor may suddenly dampen the demand 

15 for the products produced by the manufacturing plant. Thus, previously projected 
demand for products may not reflect the current and forecast demand. 

In addition, priorities originally assigned to various considerations and 
parameters on which the state function is based may have changed. The dotted line 
5508 in Figure 55 indicates a change in priorities that has arisen for the current frame 

20 of reference with respect to a prior frame of reference. For example, at the time that 
the state function was calculated, there was no knowledge of an impending road 
closure on Thursday, announced during the time interval tj - AT to ti. Knowledge of 
the road closure on Thursday may greatly increase the priority of scheduling south- 
side-resident workers for work at times when the bridge closure will least affect their 

25 commutes. 

Also, in order to prevent serious disturbances to the operation of a 
system, the system may require that previously calculated state for the system needs 
to be adhered to, regardless of whether or not a newly calculated state function would 
provide for the same state over certain intervals. For example, in Figure 55, two time 
30 windows 5510 and 5512, beginning at times T 2 and T 4 5514 and 5516, respectively, 
and ending at times T 3 and T 5 5518 and 5520, respectively, have been established, 
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within which the old value of the state function 5522 and 5524 needs to be 
maintained in a newly computed state function. Thus, historical information obtained 
in the previous state-function calculation is maintained in the current frame of 
reference. Such fence intervals may be necessary to prevent undue disruption and 
5 disturbances in system control established based on the previously calculated state 
function. 

Using the knowledge of events, knowledge of changes in 
circumstances and priorities, and the need for maintaining certain state previously 
calculated, the previously calculated and extended state function can be modified, 

10 without a de novo state function calculation, in order to produce a near-optimal state 
function for the current time window. Figure 56 shows repair of the previously 
calculated state function in order to provide a near optimal state function for a current 
time window. As shown in Figure 56, the repair of the previously calculated state 
function involves computing differences, such as difference 5602, at points in time 

1 5 within the current time window, between the previously calculated state function and 
a new, near optimal, repaired state function for the current time window. 
Computation of the differences, such as difference 5602, takes into account the 
occurrence of events in the prior interval of time tl - AT 5208, changes in 
circumstances and priorities (5508 and 5506 in Figure 55), and the need to maintain 

20 previously computed state 5522 and 5524 for particular intervals of time 5510 and 
5512. Thus, the repair-based continuous optimization technique allows for a 
somewhat continuous recalculation of state functions without full reoptimization in a 
way that minimizes disturbance and disruption of system control and prevents 
discontinuities in control. 

25 The following is a more formal description of the repair method that 

represents one embodiment of the present invention. The control optimization 
problem, solutions for which provide the control state functions that are repaired by 
the above-described repair method, are given by: 



30 



min ']Hx(t),v(t) 9 d(t))dt 
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subject to constraints: 

x l (t)-x i (t l )-)f l {x(T),v(r))dT = 0 

5 

g{x(t),v(t))>0 
veV 

f e [*„/,+ r] 

where V is a compact subset of R r and d(t) is a driving function (e.g. deterministic or 
forecasted demand). $,{f„i = \,...,n-\},an&{gjj = l,...,/w} are sufficently smooth. A 
scalar variable x n (t) > 0 is introduced such that: 

15 

X n {t)=)i{x{T), V (T),d{t))dT 

h 

*„(',) = ° 

20 The control optimization problem then becomes: 
subject to, for t e [f„f, +T) 



min x n (t, + T) 
v(,) »vi ; 



25 



*> (') - *, (<■ ) - J/ (*(*)• v{r),d(r))dT = 0 

'l 
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gi (x(t),v(t))>0 



5 where: 



v (')) fori=l,...,n-l 
' { $(x(t),v(t),d(t)) fori = n 



The solution of this problem, v (t), is called the "optimal control" and the 
10 corresponding trajectory, x*(t), is called the optimal trajectory. The above-described 
control state function (e.g. 4906 and 4916 in Figure 49) is a combination of x(t) and 
v(t), te[t t ,T). 

A previously computed control state function, with respect to the 
conventions used in Figures 52-56, is a combination of: 

v o.-ar). (/) 

x°'-^\t) 

The computed differences (e.g. 5602 in Figure 56), called "the repair," are defined to 
be a combination of: 

Sx"(t) 
Sv' : (t) 

20 Fencing indicator functions specify whether t is in an interval in which no deviations 
are allowed in the repaired control state function with respect to each control state 
variable : 



0 if t is in the "fenced" interval for x { 
1 otherwise 



25 
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0 if t is in the "fenced" interval for v k 
1 otherwise 



Assuming a finite number of fenced intervals in the interval [t l9 t t + T], then for (i = 
1,..., ri), where n is the dimension of the model, the incremental state trajectory 
5 8x[ x (t) is given by: 



Sf-(('.-A/)*)^(r)^(r) 

y=i OX j 



where the notation (f, - AT)* means that the partial derivatives are evaluated along 
1 0 the reference trajectory. 

The set {[t' M ' a 9 t'/] 9 s = l,...,Af.J is the fenced time sub-intervals for state i = 1,..., n. 

All of the events A;8(t-t e ) that happen during the time interval [? r Ar,r,] are 
approximated with a single event A} 1 (/-/,) at t x as follows: 

15 Je, (r)<fr = 'ji^J(r-/,)rfr = A? 



(r, ) = x,' 1 (f, ) - 4<~ AT) ' (t t ) - 4' provides the initial conditions for Sx? (t) 



20 Notice that v ( ''- 4r) *(r)and;c ( ' | - A7 ' ) *(/)are defined on [r, - AT, f, +T- AT]. As 
discussed above, the control state function needs to be extended to the interval 
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[J, , t x + T] . Assuming that the extension is due to a constant control, and not allowing 
fencing in the extension: 



5x'> (t x ) = 5x? (t t +T-AT)+ | 



I.+T-AT 



*=1 



« df (>,-*rr 



dd, 



{t i+ T-AT)5x){r) 
(t l+ T-*T)Sv*(T) 
{t,+T-AT)dd,{r) 



dr 



t e[t t + T-AT,t i +r],/ = l,...,n 

where <5v * (0 = v(/, +T)-v'(t, +T- AT) ' 

The dynamics on interval [t l3 f, +T] then are given by: 



10 



K(') = 



+t t ^(^ +T -^ T )*K {r)Sv) (r) 



(r) 



dr, 



/e [/„/, +r-Ar] 



i 

+r-Ar)+ J 

*, +r-Ar 



dx, 



,=1 — j 



(t,+T-AT)5x)(T) 
-(t l+ T-AT)Sv?(T) 



— 



dd, 



(t l +T-AT)Sd l (r) 



dr, 
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Expressing min x n (/, +T) as <!>(*(/, +T)) , expansion up to the second order gives: 
\(8x(*x +T)) T <M*('. +T))Sx{t x +T) 

5 

The goal of repair is to minimize the above criterion while minimizing the change 
with respect to a nominal trajectory. The combined criterion is: 

(/, + T)f <D„ (*(,, + T))fix* (/, + T) 
(/)) r 0**" (/)+±(*v* RSv* (t) 

10 where (? and R are constant positive definite matrices specified by a user that define a 
balance between satisfaction of the enterprise criterion and minimizing change. 

The constraints for the repair problem are: 
15 gj (*<''- Ar >* (t) + Sx" (/), v^* (t) + Sv" (t)) > 0J = l,...,m 
When expanded to the first order: 

g, (*<"-* r >* (/) + Sx* (/) , v^' (t) + 5v» (/)) * gj (J*-* (,), v<"->* (0) 

+^{x ( *- AT) ' (0»v ( "- Ar)> (t))Sx* (') + ff(* ( "~ Ar) * (')> v( "" Ar) * ('))*"' (r) 
= 1,...,7M 

20 

In order to satisfy the constraint inequality, given that 



gj (x^\ t )J'^\i))>0J = l,..., m 
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the second and third terms of the expanded constraint must satisfy 

5 

Meta-Control and Heirarchical Process Construction 

10 The discussion, up to this point, has focused on optimization and 

adjustment of controlling state functions generated by optimization methods at 
subsequent points in time. However, the present invention encompasses a far more 
general and powerful method for constructing complex, hierarchical computational 
systems and models, and for controlling general computation, whether simple or 

15 complex. An example of this general method for computational control has been 
presented above, when additional constraints are added to the original optimization 
problem in order to produce the meta-level, minimal-time-control-problem objective 
functional^, and the Hamiltonian for the meta-level probiem is minimized in order 
to generate a near-optimal control policy. 

20 Figure 57 illustrates a single level of meta-control of a computational 

process. In Figure 57, a first computational problem, or computational model, 5702 
appears at a lowest level along a meta-level axis 5704. The computational problem, 
or model, is discrete, and can be formally described by an iteration 
=Xk + fi( x k> u k) • As discussed above, introducing the continuous-iteration- 

25 variable parameterization of the discrete computational problem, or model, the 
original computational problem can be viewed as continuous, in order to allow for 
powerful continuous mathematical methods to be applied in order to solve the 
problem. Thus, as shown in Figure 57, the problem may be viewed as discrete, on 

one face 5706, and as continuous, ^^1 = /(jc(r), w(r)), on another face 5708. The 
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discrete-to-continuous transformation is an important enabling feature of the meta- 
control method. Using the continuous view of the problem, additional constraints 
may be added to form a new problem, such as the minimum-time control problem 
constructed in the above-described optimization problem. An important 

5 aspect of this ability to add constraints is that the newly added constraints are not 
limited or necessarily tied to the original problem. Again, in the above-described 
embodiment, the initial problem involves finding an optimal control solution for an 
optimization problem, but the added constraints relate to finding near-optimal control 
policies for steering the computation towards a time-efficient computation of the 

10 original optimization problem. Other types of constraints might be added. For 
example, one might wish to add meta-level constraints dealing with days of the week 
when optimizations should be performed, or constraints that steer the computation 
towards both time and memory-usage efficiencies. The added constraints thus 
represent a meta-level view superimposed on the original problem. The ability to add 

15 new constraints facilitates a hierarchical computational approach to solving problems, 
similar to the levels of conceptual design used in designing complex software - such 
as the well-known top-down approach to problem solving. 

In Figure 57, addition of constraints is shown as a first part of a 
transformation represented by the large arrow 5710. Notice that the entire original 

20 problem is carried forward into the new problem both within the constraints for the 
new problem and within the terms of the new objective function. Thus, the new 
problem includes the original problem, and solution of the original problem is a 
subset of the solution of the new problem. Thus, although at the meta-level, one can 
focus on new constraints and an entirely different view of the original problem, the 

25 original problem is encoded in the new problem. In a second part of the 
transformation, the Hamiltonian is used in order provide an expression for 
determining a discrete computational method for computing a meta-control policy. 
The discrete meta-control policy expression 5712 can itself be viewed as a continuous 
problem 5714 with introduction of a new continuous iteration variable o. 

30 Figure 58 illustrates multiple levels of meta-control and/or hierarchical 

construction of a complex computational process. As shown in Figure 58, repeated 
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meta-level transformations may lead from an original, discrete problem 5802 through 
an arbitrary number of meta-level problems along a meta-level axis 5810. As 
described above, each new meta-level allows for addition of new constrains and 
creation of new problems, while carrying forward the original problem and all lower 
5 meta-levels into the higher-level meta levels. As an example, consider a 
computational modeling problem of modeling a living cell, a current area of much 
interest in the scientific community. The cell may be viewed hierarchically as a large 
number of conceptual levels, beginning with a molecular level governed by quantum 
mechanics and statistical mechanics, and ascending through supermolecular 

10 complexes, macromolecular structures and compartments, organelles, genetic control, 
environmental influences, and myriad other conceptual levels, each viewed in a 
different way and as being driven by different kinds of forces and principles. 
However, each level contributes to the whole of the cell, and the effects of molecular- 
level events and conditions may profoundly influence higher conceptual levels. For 

15 example, spontaneous conversion of a particular nucleic-acid base to a different 
nucleic acid base within a chromosome may lead to an altered gene that, in turn, may 
lead to an altered gene product that is a component of a crucial organizing structure or 
temporal control mechanism that, in turn, may severely impact the viability of the 
entire cell. Moreover, various portions of the problems at each level may be naturally 

20 expressed in integer, Boolean, or other types of variables, while other portions of the 
problems may be defined in terms of continuous variables. 

The methods of the present invention are applicable to developing a 
complex computational model for a system such as a living cell or a real-time 
distributed operating system. The methods discussed above for transforming non-real 

25 variables to real-number variables and then to floating point variables provides a way 
to incorporate the natural expression of different sub-problems into an all- floating- 
point problem space to which the methods of the present invention may be applied. 
The ability to define different meta-level views of the cell-modeling problem, while 
carrying forward the lower-level views, allows for a hierarchical modeling that 

30 preserves the full complexity of lower-level views without requiring the lower-level 
views to be explicitly revisited at higher meta-levels. 
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More formally, assume that an iteration algorithm is provided as 

follows: 

where are control variables and x k are decision variables. 
5 The algorithm converges to 

lim x. = x 

k-x*> 

The specifics of the algorithm are encoded in the function f. and the parameters u k . 
In a continualization procedure, the iteration is converted to a controlled differential 
equation in terms of jc(r)and w(r)such that lim x(t) * 3c = lim x. . The variable r is a 

10 continuous iteration variable corresponding to the discrete iteration variable k. 
The differential equation is of the form: 

^)=/(*(r), W (r)) 

for all w() e U(), where U() is the set of feasible paramter values. 
The optimal control problem is formulated as follows: 

15 



subject to —£l = f(x(T),u(r)) 
w(r)et/(r) 

The functions ^ and *P are chosen by the algorithm designer in order to achieve meta 
objectives. 



Using the minimum principle of Pontryagin, and with the Hamiltonian given by 

25 

K(^(r) ) «(r), jP (r)) = ^(x(r) )M (r)) + p(r) r /(^(r),«(r)) 
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where p(x) is the costate , the necessary conditions for x*(r),w*(r),G[0,r)to be 
optimal are: 



dx'jx) 
dx 



(dH(x(x),u(x),p'(x))) 
dp(x) 



d P '{r) 
dx 



dH(x(x),u(x),p(x))) T 
dx(x) 



PK ) dx(x) 



10 



**(0) = * 0 



«*(r)eU(r) 

where the Hamiltonian has an absolute minimum as a function of u at u'(x) : 
15 H(x'(x),u{x),p(x))<H(x(x),u(x),p(x)), Vu{x)eU(x), x e[0,T) 
A convergence sequence {u s (r)} is generated by: 



{r ) = u s (x) + W(x(x),u s (x),p(x)) 



20 



such that 



limw,(7-) = «*(r), when jc(r) = x*( ), x <=[0,T) 
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The above convergence sequence generation expression is analogous to the original 
iteration algorithm 

where u s (r) is analogous to x k . 

Using the continualization procedure described above with a being a continuous 
parameterization of s, one obtains: 



10 



%^ = ^(,(r),v(r),^(r)) 

where v(r,cr) is a continulaized version of u s (t) and fF(x(r),v(r),/?(r)) is a 
continualized version of W(x(r),u s (r),/?(r)) . 



To summarize: 



15 



=/(,(,),,(*)) 



3tt(*(r), M (r),/7(r)p r 



20 



^l = IF(x(r),v(r) > p(r)) 



with /?(r) = W andjc(0) 

&(r) 



given. 
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Although the present invention has been described in terms of a 
particular embodiment, it is not intended that the invention be limited to this 
5 embodiment. Modifications within the spirit of the invention will be apparent to 
those skilled in the art. For example, an almost limitless number of different 

The foregoing description, for purposes of explanation, used specific 
nomenclature to provide a thorough understanding of the invention. However, it will 
be apparent to one skilled in the art that the specific details are not required in order 

10 to practice the invention. In other instances, well-known circuits and devices are 
shown in block diagram form in order to avoid unnecessary distraction from the 
underlying invention. Thus, the foregoing descriptions of specific embodiments of 
the present invention are presented for purposes of illustration and description; they 
are not intended to be exhaustive or to limit the invention to the precise forms 

15 disclosed, obviously many modifications and variations are possible in view of the 
above teachings. The embodiments were chosen and described in order to best 
explain the principles of the invention and its practical applications and to thereby 
enable others skilled in the art to best utilize the invention and various embodiments 
with various modifications as are suited to the particular use contemplated. It is 

20 intended that the scope of the invention be defined by the following claims and their 
equivalents: 
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